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ABSTRACT 


An  efficient  automatic  procedure  is  given  for  evaluating  the  integral  of  the  bivariate  normal 
density  function  (IBND)  over  an  arbitrary  polygon  11.  The  polygon  II,  defined  hy  N  points,  falls 
into  one  or  more  of  the  following  classes:  {3},  simple  polygons;  { S > ,  limit  elements  of  sequences 
of  uniformly  bounded  N-sided  simple  polygons  of  the  same  orientation;  (ID.  arbitrary  polygons, 
which  includes  self-intersecting  (SI)  ones,  where  {$)  C  { S)  C  (11).  It  is  not  necessary  to  specify 
the  class  beforehand.  The  method  extracts  from  n  a  set  of  N  exterior  angular  regions.  The  IBND 
is  evaluated  over  each  of  these,  and  the  results  are  properly  combined  to  yield  IBND  for  II.  In 
case  II  is  SI.  account  must  be  taken  of  the  number  of  its  "primary  circuits"  and  their  orientations. 
A  by-product  of  the  analyses  is  the  evaluation  of  a  function  A(JI>  for  which  |A|.  when  properly 
interpreted,  gives  the  area  of  II. 

Another  procedure  for  obtaining  the  same  final  results  is  described  for  completeness  which  is 
not  as  efficient.  It  treats  an  SI  polygon  by  decomposing  it  into  a  finite  set  of  S  or  S  type  elements. 
Tire  IBND  is  evaluated  over  each  of  these:  tire  results  arc  properly  summed  to  give  the  IBND  for  It. 
In  contrast  to  the  first  method,  the  smallest  class  (S),  (S).  (II)  to  which  II  belongs  must  be 
specified  for  computational  efficiency. 

The  fortran  IV  programs  for  both  procedures  are  presently  set  to  yield  approximately  d. 
or  ‘Mcdmal-digit  accuracy.  Fortran  IV  listings  of  the  programs  are  given. 
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I.  INTRODUCTION 


This  report  describes  two  automatic  and  efficient  procedures  for  evaluating  the  integral  of  the 
general  bivariate  density  function  over  an  arbitrary  polygon1  ft.  Specifically,  we  evaluate  P(fi),  i.e., 


(1) 


POI)  = 


(1  ~P2)~1/2 
2irowoz 


2  (W-Iiw)(z-Mz) 
P  OwOz 


+ 


dw  dz , 


where  ( p w,  pz)  is  the  mean  and 


potaw 


is  the  covariance  matrix  of  the  normal  random  variable  (w,  z)  with  correlation  coefficient  p,  (pi  <  1. 

Three  main  classes  of  polygonal  elements2  {$}.  { S),  {11}  are  treated  in  the  text.  The  set  of 
simple  polygons  is  denoted  by  (S)  with  the  subset  of  convex  polygons  denoted  by  {C}.  The  class 
{.$}  is  enlarged  to  include  elements  which  arc  limits  S  of  sequences  of  uniformly  bounded  N-sided 
simple  polygons  with  the  same  orientation,  (S„(N)}.  This  class  is  denoted  by  {§).  A  more  ex¬ 
tended  class  {11}  is  obtained  by  adding  self-intersecting  (SI)  polygons2  to  (S). 


Using  the  well-known  linear  transformation 


(2) 


x  * 


w  “ 


P  i-'A 

“  n 


/ 


s/l  -  p2  ,  y 


*  “  p* 


IpI  <  1 


in  (I)  results  in  a  new  integrand  which  has  circular  symmetry  about  the  origin.  In  addition,  since 
(2)  maps  straight  lines  into  straight  lines,  ft  transforms,  by  (2),  to  another  polygon  II  of  the  same 
class.  Thus  ( I )  can  be  written  as 


(3)  P(1I)  »  PUI)  -  (J  Z(x,  y)  dx  dy  t 

11 

1 A  polygon  or  polygonal  clement  will  always  mean  a  dosed  finite  broken  line,  with  it*  interior,  in  the  plane.  The 
last  segment  terminates  at  the  first  point.  Its  boundary  is  defined  by  an  ordered  set  of  N  points  in  the  plane.  How- 
ever,  wc  specify  a  polygon  by  N  ♦  1  points,  where  the  (N  ♦  1  )st  and  first  points  arc  the  same. 

2For  case  oflanguagc,  polygon  and  polygonal  dement  arc  used  interchangeably. 

JWc  say  a  polygon  is  self-intersecting  if  it  is  not  in  {S).  A  characterization  is  given  ut  Section  111. 


1 


where 


(4)  Z(x,y)  -  ^  exp(-(x2 +y2)/2]  . 

Hereafter,  we  assume  (1)  has  been  transformed  to  (3),  and  we  deal  only  with  (3),  the  integral  of 
the  bivariate  normal  density  function  (1BND)  for  II.  Also,  unless  noted  otherwise,  we  denote  an  ele¬ 
ment  of  a  particular  class  or  set  by  the  letter  in  braces  designating  that  set.  For  example,  C  refers 
to  an  element  of  (C).  Note  that  {C}  c  {S}  C  {S>  c  {II}. 

We  make  the  convention  that  if  a  simple  polygon  S  is  positively  oriented,  (PO),  i.e.,  with  its 
area  on  the  left  as  one  traverses  the  boundary  continuously,  (3)  yields  a  positive  result,  whereas  if 
S  is  negatively  oriented  (NO),  i.e.,  with  its  area  on  the  right  as  the  boundary  is  traversed  contin¬ 
uously,  (3)  yields  the  same  result  with  a  minus  sign.  If  II  is  SI,  there  can  be  both  positive  and 
negative  contributions  to  P(I1).  For  example,  in  Figure  1  below,  P(I1)  is  made  up  of  the  sum  of 
the  probabilities  over  triangles  A  and  B,  where  A  is  specified  by  the  point  set  (1,  2, 3,  1)  and  B  by 
(3,  4,  5,  6).  Thus  Pf II)  =  P(A)  +  P(B),  where  P(A)  >  0,  but  P(B)  <  0.  Clearly  then,  P(!l)  may  be 
negative  and  make  no  sense  in  terms  of  probability;  nevertheless  we  often  call  P  a  probability, 
regardless j>f  its  sign,  for  ease  of  language.  In  Figure  2,  we  have  an  example  of  an  element  S  of  {S }. 
where  P(A)  and  P(B)  are  both  positive.  These  concepts  will  be  discussed  more  fully  in  Sections  HI 
and  IV. 

in  Figure  3  we  show  an  element  of  a  sequence  (S„(  1 1 )}  of  1 1  'Sided  simple  polygons  for  which 
a  limit  element  S  is  obtained  by  letting  the  points  3-9  converge  onto  the  same  straight  line  as 
shown  in  Figure  4. 

A  main  objective  of  this  report  is  to  describe  and  discuss  two  procedures.  (P~  A)  and  (P-B). 
to  evaluate  (3).  In  cither  case,  if  H  is  in  (Sj.  P(II)  is  found  by  integrating  (4)  over  a  set  of  exterior 
angular  regions  of  U.  essentially  in  the  same  way  as  described  in  (21 .  (3)  for  convex  polygons.  The 
details  are  given  in  Section  HI.  For  background  and  completeness,  the  convex  case  is  summarised 
in  Section  II. 


S 


Figure  I .  An  SI  Polygon  of  Class  (I  I ) 


Figure  2.  A  Polygon  of  Class  (S> 


The  same  is  true  il  II  is  in  {S}.  Some  pre-processing  of  II,  which  reduces  it  to  an  element 
of  (S),  as  described  in  Section  III,  may  b-  necessary  with  (P-A)  and  desirable  with  (P-B).  It  is 
then  treated  as  indicated  in  the  previous  paragraph. 

If  II  is  self-intersecting  (SI),  procedures  (P-A)  and  (P-B)  are  quite  different.  For  (P-A),  H 
is  first  decomposed,  by  separate  sub-procedures,  into  a  set  of  disjoint  isolated  elements  in  {S}.  or 
perhaps  (S).  The  value  of  P  for  each  of  these  isolated  elements  is  computed  and  the  results  are 
summed  to  yield  P( II).  For  (P-B)  no  decomposing  procedure  is  necessary.  It  treats  11  as  if  it 
were  in  (S)  or  (S).  This  is  possible  because  it  keeps  track  of  the  number  of  “primary"  circuits,  or 
loops,  in  11.  Generally  (P-B)  is  (aster  than  (P-A),  since  there  is  no  necessary  pre-processing  to  do. 
Moreover,  in  (P-A).  for  efficiency,  one  must  specify  the  smallest  class  to  which  II  belongs: if.  by 
error,  one  specifies  a  class  to  which  II  does  not  belong,  then  P(J1)  will  be  computed  incorrectly. 
For  (P-B),  the  smallest  class  need  not  be  specified,  and  incorrectly  computed  results  are  not 
possible.  Since  (P-  A)  decomposes  II.  it  is  useful  in  analyzing  the  configuration  of  II.  Nevertheless. 
(P-B),  with  some  of  its  parts  in  common  with  (P-A).  is  the  preferred  procedure.  It  will  he  dis¬ 
cussed  in  Section  IV.  The  remainder  of  the  discussion  on  (P-A),  with  some  numerical  results,  is 
relegated  to  Appendix  A. 

\vc  emphasize  that  the  procedures  described  in  Sections  11!  and  IV  lead  to  a  computer  program 
which  is  completely  automatic  in  the  seme  that  one  can  simply  specify,  as  input,  tin?  coordinates 
used  to  define  11  in  proper  order,  the  number  N  of  such  points,  and  one  of  three  accuracies  desired 
for  P(!l).  A  by-product  of  the  program  is  a  function  A(ll),  where  |AL  properly  interpreted,  gives 
the  area  of  11. 

In  the  most  general  case  wltere  II  is  SI.  we  do  not  know  of  another  program  to  compute  P 
Kveu  if  11  is  simple  we  are  not  aware  of  an  automatic  program,  although  such  a  program  would 
liavc  many  applications  in  probability  and  statistical  studies.  Of  course,  brute  force  methods  can 
always  Ik  devised,  such  as  decomposing  It  into  triangles  and  quadrilaterals  |5.  p.  9561.  t-ven 
though  it  is  easy  to  obtain  a  program  for  decomposing  an  arbitrary  polygon  into  a  set  of  triangles, 
the  required  number  of  such  triangles  would  result  in  an  inefficient  procedure.4  For  example,  in 


4ln  Section  VI.  (see  page  47).  well  a  procedure  is  described  Combined  with  Drczncr's  algorithm,  |2.  page  18|  tl 
gives  us  an  independent  method  for  checking  out  programs.  A  listing  of  the  checkout  program  ts  given  in 
Appendix  G. 
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the  case  of  an  N-sided  polygon,  P  would  be  required,  in  general,  for  3 f N  —  2 )  angular  regions, 
whereas  (P-B)  needs  P  for  only  N  angular  regions.  The  major  time-consuming  portion  of  the 
program  for  computing  P(  FI)  is  the  evaluation  of  P  for  each  angular  region  needed. 

The  computer  program  and  its  flow  charts  are  discussed  in  Section  V.  The  Fortran  IV 
CDC-67005  program  listings  are  given  in  Appendix  F.  Section  VI  contains  numerical  results  for 
polygonal  elements  displayed  in  Figures  30-58.  These  figures  demonstrate  the  robustness  of 
the  program. 


II.  NORMAL  PROBABILITY  OVER  CONVEX  POLYGONS  (SUMMARY) 


The  numerical  method  and  computing  program  for  evaluating  P(C),  where  C  denotes  a  convex 
polygon,  are  described  in  detail  in  an  earlier  report,  [2] .  We  proceed  to  summarize  the  main  ideas 
since  most  of  them  carry  over  to  polygons  in  (S) . 

We  use  the  notion  of  an  angular  region,  i.e.,  the  semi-infinite  region  bounded  by  two  inter¬ 
secting  straight  lines.  There  are  four  such  angular  regions  associated  with  two  intersecting  straight 
lines  and  one  must  keep  in  mind  which  one  is  of  interest.  Let  a  denote  the  angular  region  of 
interest.  As  shown  in  Figure  5.  it  can  be  specified  by  the  parameters  R,  0,.  0:,  Lines  I  and  2 
form  the  boundaries  of  a.  The  quantity  R  denotes  the  distance  from  the  origin  to  the  vertex  of 
a  at  V.  When  necessary  we  shall  denote  a  by  a(R.  0,,  0,).  Our  objective  is  to  give  an  efficient 
procedure  to  evaluate  P{ a). 

Because  of  the  circular  symmetry  in  the  integrand  of  (3).  it  is  convenient  to  perform  a  rota¬ 
tion  of  axes  through  the  angle  r,  such  that  the  line  L  and  the  new  positive  x-axis  coincide.  The 
rotation  for  Figure  5  is  shown  in  Figure  6.  Hereafter,  we  assume  such  a  rotation  has  been  carried 
out  and  call  the  new  axes  x  and  y  again. 


5  the  CtX'-o?00  u  i  taige-scale  binary  computer  which  doc*  10*  arithmetic  operations  per  second  on  4&-biI  mammas. 
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y* 


Figure  6.  Showing  Angular  Region  a 
Relative  to  Rotated  Axes 


Introducing  polar  coordinates  (r,  0)  centered  at  V  =  (R,  0),  we  have 

(5)  x  =  R  +  r  cos  0  ,  y  =  r  sin  0  ,  -ir  <  0  <  ir . 

The  variable  0  will  be  measured,  from  the  x-axis  (x  >  R)  about  (R,  0),  as  positive  in  the  counter¬ 
clockwise  direction.  Using  (5),  (3)  becomes,  with  a  in  place  of  FI, 

(6)  P(a)  =  '2n  j  f  exP  [~ 1/2(R2  +  2rR  cos0  +  r2)]  rdr  d0  . 

d\  0 

An  integration  by  parts  in  (6),  on  the  integral  in  r,  yields 


(7) 
where 

(8) 


e-W2e-rRcos^  rdr  =  1  -  2u  erfc(u)/z(u) , 


_  R 


u  = 


VT 


COS0,  z(u) 


*4=  exp  (— u2) , 
V* 


erfc(u)  s 


z(t)  dt . 


Using  (7)  in  (6),  and  carrying  out  the  obvious  part  of  the  0-integration  yields 


(9) 


P(a)  =  e-Rl/2 


u(erfc(u)/z(u)]  d 0 


If  R  =  0,  then  P(a)  =  A0/2ir  as  required,  where  A0  =  02  -  0, . 

For  exterior  angular  regions  of  polygons  we  can  require  that  -rr  <  A0  <  ir.  For  PO  (NO) 
convex  polygons,  02  will  always  precede  0*  in  the  counterclockwise  (clockwise)  direction,  so  that 
0  <  A0  <  ir  (— 7r  <  A0  <  0).  Hence  from  ((>},  P(o)  >  (<)0  for  PO  (NO)  convex  polygons. 
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We  resolve  the  difficulty  of  evaluating  the  integral  in  (9)  by  using  a  minimax  polynomial  fit  to 
erfc(u)/z(u)  for  0  <  u  <  c(5).  That  is,  for  a  given  6  >  0,  a  set  of  real  numbers,  {ak},  and  a  least 
positive  integer  K  are  found  such  that 


(10) 


K 

erfc(u)  -  z(u)  ^  akuk 
o 


0  <  u  <  c(5). 


It  is  shown  in  [2,  page  6]  that  if  ( 10)  holds  then 
-R2/2  rO 2  f 


(11) 


R2/2  CV2  \  ^  ^ 

—  j  ju(erfc(u)/z(u)]  -  akuk+1 


d0 


<  6/y/W. 


Thus  the  constant  c  is  chosen,  once  5  is  specified,  so  that  the  probability  over  the  semi-infinite 
region  to  the  right  of  the  line  x  =  \/2  c  is  equal  to  6/s/jT,  i.e., 

(12)  \  erfc(c)  =  e  ^  6/>/F 


The  coefficients  ak  as  well  as  K  and  c(6)  are  given  in  [21  for  4  different  values  of  e  corresponding 
to  desired  accuracies  of '  6,  9,  and  12  decimal  digits  in  P(o).  They  are  also  included  in  Appendix  E 
of  this  report. 

Recurrence  relations  allow  us  now  to  carry  out  the  numerical  integration  of  the  integral  in  (9). 
We  have,  using  (10),  that  P(a)  is  given  within  ie  by 


(13) 
where 

(14) 

with 

(15) 

and 

(16) 


P(a)  = 


e-RJ/2 


AO 


it/  V'  , 

2  ~  Li  akJk*l 


,  10,1  I02i 


r 


< 


r0  j 

Jk  3  {W/Vzf  cosk0dO 

Jk*l  “  k~+-i  “Ml)  +  k  (3)  4-1  ]> 


R 


1) 

h,  =>  sin  0, ,  i  =  1.2 


6t  n  cos  &i 
R2  =  (x2  +  y2) ,  with  vertex  of  a  at  (x.  y) , 


J0  =  A0  =  P2  -  0,  .  J,  *  h2  -  h,  . 
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The  constraints  |0j  I  <  it  12,  |02 1  <  ff/2  in  (13)  are  required  since  (10)  only  holds  for  0,  which 
requires  cos  6  >  0,  since  R  >  0.  For  arguments  outside  the  range  \6\  <  it/ 2,  we  make  use  of  the 
relation 


(17)  P[a(R,0,fl)]  =  j  erfc  (^sinfl)  -  P[a(R,O,ir~0)],  \6\  <  it , 

where  a(R,  0,  6)  denotes  the  angular  region  with  its  vertex  at  R,  0j  =  0,  02  =  #  (see  Figure  6). 
The  various  situations  in  which  (1 7)  is  needed  are  shown  in  Figure  5  of  [2] . 

The  program  resulting  from  (13),  (14),  (17)  is  very  efficient  and  takes  advantage  of  situations 
where  reduced  computing  effort  is  possible,  namely  when  R  is  sufficiently  large  or  small.  For 
example,  ( 1.8)  is  used  when  R2/2  <  a2,  with  G  set  to  zero  when  R2/2  <  oq. 


(18) 


S  f  -  G  *  f 


LaVF 


(h2  ~  hl>  "  2»  (82h2”glh!) 


See  Section  V  (page  28)  for  added  comments.  Details  are  given  in  [2,  pp.  6-8] . 

Letting  C(v(,  v2 . vN)  denote  a  convex  polygon  with  N  vertices,  vlt  v2 . vN,  where 

(Vj)  s  (xj,  yj),  we  have 

N 

(19)  P(C)  =  1  -  £  P(flj)- 

1 


litis  equation  is  fundamental  to  computing  P(C)  by  the  use  of  probabilities  over  appropriate 
angular  regions  at.  i  a  l,  ....  N.  These  angular  regions,  quite  simply,  turn  out  to  be  the  exterior 
angles  of  C  as  shown  in  Figure  7  for  N  a  6. 6 


Figure  7.  Convex  Polygon,  N  «  6.  Shows 
Angular  Regions  for  (19). 


6ht  (2] .  (3)  we  designated  die  angular  regions  for  P(C)  in  a  slightly  different  but  equivalent  way. 
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Ciearly  (19)  is  true  since  2^  P(flj)  =  1  -  P(C).  We  also  note  that  the  vertices  are  specified  in 
counterclockwise  order  so  that  the  area  of  C  is  on  the  left  as  one  travels  along  C  continuously  from 
0\)  to(vk,.,),  k  =  1,  N,  (vN+1)  s(V[). 

When  A 6  is  very  near  0  or  ir,  the  possibility  of  a  catastrophic  error  due  to  round-off  must  be 
dealt  with.  As  an  example  of  this  singular  situation,  suppose  we  are  considering  a  polygon  where 
one  of  the  angular  legions,  say  a,  as  shown  by  the  solid  lines  in  Figure  8,  subtends  an  angle  of 
nearly  n  radians  with  the  sides  of  a  at  large  perpendicular  distances  from  the  origin,  so  that 
P (a)  ~  1 .0.  But  suppose  that  by  rounding  error  line  (7)  is  actually  given  by  the  computer  as 
line  (5)  so  that  the  angular  region  a  subtends  an  angle  0,  near  (~ir)  radians.  The  program  would 
then  yield  a  value  P(a)  near  zero.  Moreover  it  would  be  negative  since  0  is  measured,  from  (3)  to 
(T),  clockwise  rather  than  counterclockwise  as  required.  The  reader  should  note,  as  stated  earlier, 
that  for  a  positively  oriented  convex  polygon  all  the  angular  regions  Oj  i  =  1,  ...,  N  are  positive, 
i.e.,  rotating  from  6X  to  02  is  always  in  the  counterclockwise  direction  so  that  each  A 0;  is  non¬ 
negative  and  no  larger  than  rr. 

The  program  is  alerted  to  a  singular  situation,  such  as  the  above  example,  whenever 

(20)  sin(02  -0t)  <  0, 

for  any  angular  region  a,,  ....  aN.  If  thic  inequality  holds,  and  it  can  only  hold  through  rounding 
error,  a  second  inequality  is  tested,  name!;. 

(21)  cos'02  -0,)  <  C. 

If  (20)  and  (21)  arc  both  satisfied,  we  resolve  the  difficulty  by  setting  P(C)  -  0,  since  if  an  angular 
region  has  A0  =  ?r  for  a  convex  polygon  all  vertices  must  be  on  the  same  line.  If  (20)  holds  but  (21) 
does  not,  we  set  P(tr)  =  0,  since  AO  ~  IQ-14.  However,  in  this  case,  a  remote  possibility  can 
arise  for  which  P(a)  =  0  mav  be  incorrect.  In  particu'ir,  when  g,  and  g,  are  both  negative  and  R  is 
verv  large,  say  109,  P(a)  may  not  be  near  zero  even  though  AO  ~  10-'4.  Under  these  very  unlikely 
circumstances,  we  cannot  determine  P(a)  within  the  single  precision  capabilities  of  theCIX?*fi?00, 
because  of  the  inadequate  precision  in  A 0. 


Figure  8.  A  Singular  Case  Situation 
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In  the  earlier  report  on  convex  polygons,  [2],  we  remarked  that  a  program  was  available  for 
arbitrary  polygons  which  made  fundamental  use  of  the  probability  program  for  convex  polygons. 
Those  ideas  are  briefly  summarized  in  Appendix  C.  Since  the  completions  of  [2] ,  [3]  later  studies 
revealed  that  a  significant  increase  in  efficiency  could  be  made  by  dealing  directly  with  simple 
polygons,  rather  than  decomposing  them  into  sets  of  disjoint  convex  polygons,  as  done  originally. 
These  results  are  the  main  topic  of  the  next  section. 


III.  NORMAL  PROBABILITY  OVER  S  AND  S  POLYGONS 

In  this  section,  we  describe  our  method,  based  on  computing  normal  probabilities  over  exterior 
angular  regions,  for  evaluating  P  for  elements  in  {S>.  The  same  method,  by  continuity  arguments, 
can  be  used  to  find  P  for  elements  in  {S},  provided  certain  precautions  are  taken,  which  will  be 
discussed  later.  Just  as  for  the  class  (C>  (convex  polygons),  we  shall  show  that  P  for  an  N-sided 
simple  polygon  requires  the  integration  of  (4)  over  its  N  exterior  angular  regions.  Taking  the 
precautions  mentioned  above  into  account,  no  more  than  N  integrations  are  also  needed  to  compute 
P(S),  where  S  is  specified  by  N  points. 

It  is  important  to  keep  in  mind  that  the  angular  measure  AO  =  02  ~  0y  for  an  exterior  angular 
region  a  of  a  polygon,  satisfies  -v  <  A0  <  ir.7  Also  one  can  see  from  (6)  that  P (a)  takes  the  same 
sign  as  AO,  where  0  is  taken  positive  (negative)  measured  from  the  x-axis,  x  >  R,  about  (R,  0)  in 
the  counterclockwise  (clockwise)  direction,  (see  (5)). 

We  say  S  is  positively  oriented  (PO)  if  its  vertices  (vj),  j  =  1 ,  2 . N  are  ordered  such  that  the 

interior  of  S  is  on  the  left  as  the  boundary  is  traversed  continuously  in  the  direction  of  increasing  j. 
If,  on  the  other  hand,  the  interior  of  S  is  on  the  right  as  the  boundary  is  traversed  in  the  way  just 
described,  then  we  say  S  is  negatively  oriented  (NO).  For  example,  the  polygons  A  and  B,  making 
up  II.  in  Figure  l  are  1*0  and  NO,  respectively. 

One  way  to  determine  the  orientation  of  S  is  by  the  sign  of  the  expression  for  A(S). 

1  n 

ACS)  a  j  £  Xj(yj4,  -yH).  y0  a  yN,  yN>,  s  yt, 

J-i 

where  (Xj.  yj)  denotes  the  coordinates  of  the  j,h  vertex  (v})  of  S.  In  fact,  it  can  be  shown  that  S  is 
1*0  if  and  only  if  A(S)  >  0,  and  NO  if  and  only  if  A(S)  <  0,  (see  Appendix  D).  The  area  of  S  is 
given  by  I A  I.  The  expression  in  (22)  yields  an  efficient  procedure  for  computing  A (S).  The  com¬ 
puter  program  for  this  computation,  SMP*7,  is  listed  in  Appendix  F.  The  derivation  and  orientation 
properties  of  (22)  are  given  in  Appendix  D.  Also  the  expression  for  A  in  the  wz-plane  is  given  there, 
(See  page  1 ). 

Consistent  with  earlier  remarks  that  P(S)  is  taken  positive  if  S  is  1*0,  we  have 
flHS)  >  0.  if  A  >  0. 

|l*(S)<0,  if  A  <  0.  (P(S)  =  0,  if  A  a  0).  (See  page  26). 

7 For  compuutiom,  a  strict  inequality  on  die  left  is  used.  More  on  litis  later  in  the  report. 
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Moreover,  since  A  is  a  continuous  function  of  the  vertex  coordinates  of  S,  (22)  must  also  hold  for 
all  S  in  {S}.  Thus  S  may  be  replaced  by  S  in  (23). 

Consider  first  that  S  is  PO;  then  we  shall  show  that 

N 

(24)  P(S)  =  1  P(fli). 

i=l 

Here  av  as  before,  denotes  the  exterior  angular  region  at  the  i,h  vertex  of  S,  which  is  formed,  as 
always,  by  extending  the  sides  (i  -  1,  i)  and  (i,  i  +  1)  as  shown  in  Figures  7  and  9,  where  for  i  =  1, 
(0,  1)  =  (N,  1)  and  for  i  =  N,  (N,  N  +  1)  =  (N,  1).  A  glance  shows  that  (24)  is  the  same  as  the 
expression  for  P(C)  given  by  (19).  In  (19)  each  P(aj)  is  positive.  In  (24)  this  will  not  be  the  case  if 
the  interior  angle  at  (vt)  of  S  exceeds  tt  radians.  For  example,  in  Figure  9,  the  interior  angle  at  (3) 
exceeds  ir,  so  that  a3  is  measured  in  the  clockwise  direction  rather  than  counterclockwise.  Hence 
Ad  <  0  for  a3  and  P(a3)  <  0.  Note  that  P(flj)  >  0  for  i  =  1,2, 4. 


Figure  9.  Polygon  S.  Shows  Angular 
Regions^,  i=  1, ....  4. 


Since  the  sign  of  P(a)  is  determined  by  the  sign  of  A0t  it  may  also  be  fixed  by  the  sign  of 
sin  AO.  Thus,  we  also  have 


P(«)  >  0  if  sin  AO  >  0, 

(25) 

P(a)  <  0  if  sin  AO  <0.  AO  &  03  -  0, ,  IA0I  <  *. 

As  a  second  example,  in  Figure  3,  on  page  3,  S  is  PO  with  P (a { )  >0  for  i  “  1,  2.  3,  5.  7, 9,  10,  II, 
and  P(a,)  <  0  for  i  =  4,  6,  8.  The  ambiguous  case,  sin  A0  »  0,  with  \A0\  =  rr,  is  connected  with  the 
precautions  mentioned  earlier  for  S,  and  will  be  discussed  later  in  this  section. 


Let  (S+)  denote  S  when  PO,  and  let  (S-)  denote  the  same  configuration  when  NO.  If  (24) 
holds  for  (S+),  then 

N 

(26)  P(S-)  «  -1  P(fl,). 


where  uNtI.(  from  (24)  and  a,  from  (26),  with  (N  +  l)  =  ( I),  are  vertical  angles  with  their  measures 
of  opposite  sign. 
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Indeed,  since  P(S-)  =  -P(S+)  =  -I  +  2  P(aj),  and  2  P(flj)  for  (S-)  has  the  same  absolute  value, 
but  opposite  sign  as  the  corresponding  quantity  for  (S+),  (26)  follows  directly.  By  continuity 
arguments,  (26)  also  holds  for  NO  elements  in  (S)  provided  certain  precautions  are  taken. 

The  truth  of  (19)  for  convex  polygons  is  obvious.  In  the  case  of  (24)  we  give  a  heuristic 
argument  for  its  validity,  which  can  be  made  rigorous. 

The  argument  is  inductive.  Certainly  (24)  holds  for  N  =  3.  Some  insight  is  gained  by  con¬ 
sidering  N  =  4  with  S  not  convex,  as  in  Figure  9.  We  see  that  1  -  P(S)  is  obtained  by  considering 
£j  P(flj),  where  Pfuj),  which  is  negative,  compensates  exactly  for  the  excessive  positive  contribu¬ 
tion  from  P(flj).  Thus  (24)  holds  for  N  =  4. 

Now  assume  (24)  is  true  for  N  =  J  -  1 ,  J  >  4.  We  shall  show  that  (24)  holds  also  for  N  =  J. 
We  look  at  the  special  case  J  =  8,  with  Figure  10,  since  the  essentials  of  a  rigorous  proof  are  con¬ 
tained  in  the  arguments  for  this  case. 

First  a  diagonal  L  is  drawn  from  vertex  (3)  to  vertex  (7)  which  remains  inside  S.  Such  a  line 
can  always  be  found  for  any  simple  polygon;  a  proof  of  this  fact  is  given  in  Appendix  B.  This  line 
separates  S  into  two  simple  polygons  with  the  same  orientation  as  S.  Call  the  separated  polygons  D 
and  E.  Each  has  fewer  than  J  vertices.  From  Figure  10,  D  is  defined  by  the  vertices,  of  S,  ( 1),  (2), 
(3),  (7),  (8),  and  E  by  the  vertices  (3),  (4),  (5),  (6),  (7).  Note  S,  D,  and  E  are  all  PO.  By  assumption. 
(24)  holds  for  both  D  and  E,  where  we  have  ai(  i  ~  1,2,  ....  8,  the  angular  regions  of  S;£/j  and 
j  =  1.  2, 3, 4, 5,  the  angular  regions  of  D  and  E.  Hence 


II 


(27) 


P(D)  +  P(E)  =  P(S) , 


P(D)  =  1  P(^),  P(E)  =  1  P(^), 


i=i 


i=l 


where 

(28)  i/j  —  tf  j ,  £/*2  —  a2,  ^5  =  flg)  &2  ~  a4>  ^3  =  a5>  *^4  ~  a6‘ 

We  also  have  ^ ^ ,  and  as  shown  in  Figure  10.  Moreover,  with  P(a3)  <  0  and  P(a7)  > 0, 
clearly 

(29)  P(a3)  =  P(£/»3)  -  [P(L)  -  P(^)],  P(a7)  =  P(*f5)  +  P(94)  -  [1  -  P(L)], 

where  P(L)  denotes  the  (positive)  probability  over  the  half-plane  below  the  extended  line  L. 
From  (27), 

J  5  5 

(30)  P(D)  +  P(E)  =  2  -  £  P(£>j)  -  £  P(tfj). 

M  i«l 


By  using  (28)  and  (29)  in  (30)  wo  have 

(31)  P(S)  =  P(D)  +  P(E)  =  2  -  P(at)  -  P(a2)  -  (P(a3)  +  P(L)  -  P(tf,)] 

-(P(a7)  -  P(rf5)  +  1  -  P(L))  -  P(a8)  -  P(<f,)  -  P(a4)  -  P(a5) 
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-Pt«6>  -  P(<fj)  =1-^  P(ai)- 

This  completes  the  argument  based  on  Figure  10.  In  order  to  make  the  proof  rigorous,  it  is  necessary 
to  consider  all  the  basically  different  possibilities  for  the  measures  of  the  angular  regions  at  the  two 
vertices  on  the  diagonal  L.  In  Figure  10,  the  interior  angle  at  (3)  was  greater  than  v  radians  and  the 
one  at  (7)  less  than  tr  radians.  The  three  other  possibilities  were  checked.  The  arguments  in  these 
cases  required  nothing  new  and  are  not  included. 

In  discussing  the  sign  associated  with  P(a),  see  (25),  the  ambiguous  case  arose  where  sin  A0  a  0, 
I  AO  |  =  ir.  The  precautions  we  mentioned  at  the  beginning  of  this  section,  in  order  to  extend  our 
results  from  the  class  {S>  to  the  class  (S),  are  necessary  because  of  this  ambiguous  case  for  P(a). 
To  specify  those  precautions  explicitly,  we  introduce  several  definitions. 

Wc  say  an  angular  region  a  is  singular  (SAR)  if  (a)  or  (b)  holds: 

(a)  Three  consecutive  points  specifying  a,  say  (k  -  1),  (k),  (k  +  1)  are  colincar,  i.c.,  are  on 
the  same  tine  L,  with  |A0|  =  ».  Such  an  angular  region  is  called  a  t -angular  region, 
(PAR). 

(b)  Two  successive  points  of  a,  say  (k  -  1),  (k),  coincide.  Wc  call  them  successive  duplicate 
points,  (SDP). 
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Examples  of  (a)  occur  in  Figure  4  for  k  =  4,  5,  6,  7,  8.  Examples  of  (b)  are  shown  in  Figure  31, 
where  ( 1 ),  (2)  and  (3),  (4)  and  ( 1 5),  ( 1 6)  are  SDP. 

The  problem  jjf  SDP  is  easily  treated  by  discarding  one  of  the  points  from  the  xy-array  speci¬ 
fying  the  polygon  S.8  Clearly  dropping  such  points  does  not  affect  the  value  of  P(S). 

The  PAR,  say  a,  which  is  our  prime  concern  here,  has  the  property  that  its  angular  measure  Ad 
can  take  either  a  plus  or  minus  value  of  n  and  P(a)  given,  with  a  proper  choice  of  signs,  by  ±1/2  erfc 
(±t/\/5),  where  t  denotes  the  normal  distance  from  the  extended  line  L  to  the  origin.  Once  the  sign 
of  A8  is  chosen  “correctly.”  and  this  is  not  trivial  to  do,  the  signs  of  P(a)  and  of  the  argument  t  are 
determined.  In  addition,  from  a  computational  viewpoint,  the  correct  choice  of  signs  for  a  PAR 
must  take  into  account  the  fact  that  the  4-quadrant  arctangent  subroutine  returns  A 6  =  n  which 
may  be  wrong.  If  the  correct  choice  is  not  made,  then  computationally  such  an  element  of  {S}  is 
considered  to  be  SI.  The  reasoning  for  this  must  be  postponed  until  we  give  a  characterization  of 
SI  polygons  in  the  next  section. 

As  a  consequence  of  SAR(s),  our  main  program  package  contains  two  subroutines  for  evaluating 
P(S)  which  treat  SAR(s)  in  different  ways.  One,  VALR-7,  cannot  evaluate  P(S)  directly  if  SAR(s) 
exist  in  S.  It  resolves  the  problem  by  pre-processing  S  with  another  subroutine  SORT  III  which 
eliminates  SAR(s).  This  is  permissible  since  the  deletion  of  such  regions  does  not  affect  the  value 
of  P(S).  For  example,  SORT  III  in  processing  S  of  Figure  4,  would  delete  points  (4),  (5),  (6),  (7), 
(8)  which  clearly  would  not  change  the  value  of  P(S).  The  deletions  by  SORT  III  reduce  S  to  a 
simple  polygon  which  can  then  be  treated  by  VALR-7  to  compute  P(S). 

The  other  subroutine  for  evaluating  P(S).  VALR-2,  based  on  (P-B),  is  the  more  versatile 
routine.  In  Sections  IV  and  V.  we  shall  show  that  it  can  find  P(ri)  for  any  II  in  {II}.  It  does  not 
need  to  pre-process  II.  and  yet  it  handles  PAR(s)  and  SDP  so  that  P(I1)  is  always  computed  cor¬ 
rectly.  Thus,  it  has  no  problem  with  singular  situations  mentioned  on  page  8.  Like  VaLR-7,  it 
integrates  (4)  over  N-angular  regions  to  determine  P(II),  where  II  is  specified  by  N  points. 

The  subroutine  VALR-7  is  included  in  the  program  package,  because  it  can  be,  on  the  average, 
slightly  faster  than  VALR-2  for  simple  polygons,  where,  of  course,  most  applications  occur.  It  uses 
A(SI,  ( 22),  to  decide  whether  (24)  or  (26)  is  needed.  The  reason  for  its  greater  speed  will  be  given 
in  Section  V,  page  38, 

hi  the  next  section,  the  discussion  is  extended  to  arbitrary  polygons. 


IV.  NORMAL  PROBABILITY  OVER  ARBITRARY  POLYGONS 

In  this  section  we  show  how  to  evaluate  (3)  for  SI  polygons.  By  our  earlier  remarks,  these 
elements  belong  to  { II)  but  not  to  {S}. 


8Wc  limit  our  discussion  here  to  dements  in  (S),  but  certainly  SAR(s)  can  also  occur  with  dements  of  {11}  not 
in  {S}. 


13 


Recall  that  (P- A),  to  be  discussed  in  Appendix  A,  is  based  on  decomposing  an  SI  polygon  II 
into  a  set  of  disjoint  elements9 lo  in  {S }  or  {§}.  In  addition,  if  care  is  not  taken,  and  a  class  is  speci¬ 
fied  to  which  II  does  not  belong,  then  a  wrong  value  of  P  will  result;  moreover  for  computational 
efficiency,  it  is  necessary  to  specify  the  smallest  class  to  which  n  belongs. 

For  (P  —  B),  on  the  other  hand,  we  shall  show  it  is  not  necessary  to  specify  the  class  to  which  n 
belongs,  and  that  P  for  any  II  in  {11}  is  evaluated  by  computing  P  over  its  N  exterior  angular 
regions,  where  N,  as  usual,  denotes  the  number  of  points  specifying  II. 

We  first  characterize  SI  polygons.  Then  we  describe  in  detail  and  establish  procedure  (P-B), 
(see  page  3)  for  evaluating  P(II). 

One  way  to  note  an  SI  polygon  is  to  show  it  is  not  in  {S}.  Our  classification  procedure  is 
straightforward,  and  it  is  easy  to  conclude  from  it,  in  principle,  if  a  polygon  is  SI.  Before  supporting 
these  statements,  some  additional  definitions  and  notation  are  given. 

The  j,h  node  (j)  associates  the  integer  j  with  the  jth  xy-point,  (xj,  yj)  of  the  ordered  point 
set  V  which  defines  H  The  set  V  is  denoted  by  the  set  of  nodes  ( 1 ,  2, . . . ,  j  -  1 ,  j,  j  +  1 ,  . . . ,  N  + 1 ). 
Let  the  jth  edge  of  II,  denoted  by  j,  mean  the  directed  line  segment  of  II  originating  at  (j  -  1)  and 
terminating  at  (j),  so  that  j  terminates  and  j  +  1  (N  +  1  =  T)  begins  at  (j).  We  say  j  and  j  +  1  are 
associated  with  the  j,h  node.  Of  course,  more  than  one  node  can  exist  at  the  same  xy-point,  in 
which  case  that  point  is  called  a  multiple  node  (MN).  If  only  one  node  occurs,  it  is  called  a  simple 
node  (SN).  For  example,  in  Figure  13,  page  16,  the  first  node  has  two  edges  I,  2  associated  with  it. 
Nodes  (4)  and  (7)  at  that  same  xy-point  have  the  edges  4,  5  and  7,  8  associated  with  them,  respec¬ 
tively.  We  identify  a  particular  MN  by  MN(i),  where  (i)  refers  to  the  first  node  met  at  that  xy-point. 
In  the  example  above  of  Figure  1 3,  we  would  refer  to  MN(  1 ).  A  vertex  ( j )  of  II  is  a  point  of  V  such 
that  J  and  j  +  1  have  different  slopes.  We  define  a  path  (j,  k )  of  II  as  a  line  made  up  of  consecutive 
edges  j,  j  +  1 . k,  j  <k. 

In  order  to  characterize  an  SI  polygon,  we  use  an  oe-numbertng  scheme,  sometimes  simply 
designated  as  an  adoption,  for  specifying  II.  By  this  scheme,  all  vertices,  points  where  two  segments 
cross,  and  initial  and  terminal  points  of  overlapping  segments,  am  numbered  in  their  natural  order 
as  II  is  traced;  i.e„  starting  from  a  point  (1),  each  time  such  a  situation  is  met  it  is  numbered,  in 
sequential  order,  until  II  is  completely  traced.  Some  polygons  numbered  under  the  adoption  are 
shown  in  Figures  13,  IS,  18,  19,  20,  24,  25,  and  one  that  is  not  is  given  in  Figure  26. 10  We  will 
give  a  brief  further  discussion  on  Figures  25  and  26  at  the  end  of  this  section,  (Page  24). 

To  establish  whether  II  is  SI,  we  use  the  fact  that  II  cannot  be  a  limit  element  of  (S„(N)),  a 
sequence  of  uniformly  bounded  N-sided  simple  polygons  with  the  same  orientation,  if  the  path 
formed  by  two  of  its  edges  associated  with  a  nodal  point  at  an  MN  penetrates  another  such  path  at 
the  same  MN.  By  penetrate,  we  mean  pass  through  rather  than  just  meet.  Such  a  situation  is  shown 
in  Figure  11.  If  two  paths  (j,  j  +  1 1 .  (k,  k  +  1 1  just  meet  at  an  MN,  as  shown  in  Figure  12,  II  may 
be  reached  by  a  sequence  (Sn(N)},  as  in  Figure  2,  page  2. 


9Two  polygonal  dements  of  11  are  disjoint  if  they  have  no  more  than  one  node  in  common.  A  node  is  defined  in 
the  fifth  paragraph  on  this  page.  A  set  of  elements  is  disjoint  if  its  elements  are  pairwise  disjoint. 

loOfton,  for  computations,  die  number  of  nodes  under  the  a-option  can  be  reduced. 
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Using  these  notions,  it  is  easy  to  see  that  the  polygon  of  Figure  1  is  SI,  whereas  the  one  in 
Figure  2  is  in  {S}.  There  are,  however,  configurations  with  overlapping  edges,  such  as  appear  in 
Figures  15,  18,  and  19,  (page  17),  which  need  a  more  precise  description  for  SI  polygons. 

With  this  in  mind  and  with  MN(j)  denoting  the  MN  at  (j),  where  (j)  is  the  first  numbered 
node  at  that  point,  let  J(j,  5^  represent  a  disk  centered  at  (j)  with  radius  5j,  where  6j  is  chosen  so 
small  that  J  (j,  5p  intersects  only  those  edges  of  n  which  originate  or  terminate  at  MN(j).  Often, 
we  simply  refer  to  such  a  disk  as  a  J~disk. 

The  approach  now  is  to  make  a  T-construction,  i.e„  to  construct  a  closed  polygonal  path  T, 
which  is  “close”  to  11,  and  which  is  then  used  to  classify  II,  as  to  type,  S  or  SI,  So,  consider  two 
successive  nodes  (k  -  1 ),  (k)11  with  2  <  k  <  N  +  1,  and  (N  +  I )  s  ( 1 ).  Let  Tk  denote  a  segment  of  T 
which  will  be  taken  “close"  to  edge  k.  This  segment  is  constructed  as  follows:  Begin  with  k  =  2. 

(A)  Tk  =  k,  if  (k  -  )  )and  (k)  are  SN  points,  (N  +  1)  s  (l). 

(B)  Tk  =  Bk.  if (k  -  1)  isan  SNand (k)isat  MN(j).  1  <j  <k.  where  Bk  emanates  from  (k  -  1) 
and  terminates  at  a  point  tk  in  J(j,  5j).  If  k  *  N  +  I ,  then  we  require  j  =  l .  T ^ + j  =  T,  = 
Uj  =  T  so  that  T  is  closed,  i.e.,  T ,  terminates  and  T2  originates  at  ( 1 ). 

(0)  Tk  =  Dk,  if  (k-  1)  is  at  MN(i)  and  (k)isun  SN,  1  <  i  <  k,  where  l)k  emanates  from 
tk_j .  a  point  in  J(i,  6t),  and  terminates  at  (k).  If  i  =  1  and  k  =  2,  then  t,  =  (1)  so  that  T 
will  be  closed. 

(D)  Tk  =  Lk.  if  (k-  1)  is  at  MN(i)  and  (k)  is  at  MN(j),  i  *  j.  I  <i<k.  I  <  j  < k.  where  l  k 
emanates  from  tk_,  in  J(i.  6()  and  terminates  at  tk  in  J  (j.  5p-  If  k  =  N  +  1.  then  j  *  1. 
and  for  T  to  bo  closed,  TN+1  a  T,  =  Lj  with  tN+1  =  t1  =  (1).  Also,  for  the  same  reason, 
if  k  3  2,  i »  1  und  t(  -Cl). 


11  Wiih  no  loss  in  generality  (k  “  1 )  and  (k)  arc  assumed  not  to  be  SDP. 
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Repeat  the  procedure  with  k  =  3,  4,  ....  N  +  1  to  obtain  T.  Clearly  by  choosing  the  5j  sufficiently 
small  T  can  be  obtained  arbitrarily  close  to  FI.  Now  if  the  tk  can  be  chosen,  for  any  5j  so  that  T  is 
simple,  then  by  choosing  a  sequence  of  5;  approaching  zero,  for  each  5{,  a  sequence  of  T’scan  be 
constructed  which  make  up  {Sn(N)}  converging  to  IT.  Hence  in  this  case  n  is  in  (S } .  If  this  cannot 
be  done,  i.e.,  if  in  some  J-disk  paths  of  T  must  cross  then  n  is  SI  since  it  cannot  be  obtained  as  the 
limit  of  a  sequence  (Sn(N)}.  If  an  intersection  takes  place  in  J(j,  5j)  between  paths  [Tk,  Tk+1]  and 
[Tk+m-  Tk+m  +  i  ]  we  say  II  has  a  SI  point  at  (k,  k  +  m). 

Clearly  from  this  characterization  of  S  and  SI  elements,  by  T,  the  polygon  in  Figure  1  is  SI 
and  the  one  in  Figure  2  is  in  {S}.  A  more  complex  example,  given  in  Figure  13,  shows  it  is  neces¬ 
sary  to  consider  all  of  the  nodes  at  an  MN.  Accordingly,  it  is  not  determined  that  II  of  Figure  13 
is  SI  until  J ( 1 ,  8[ )  is  entered  for  the  sixth  time,  as  shown  in  Figure  14,  where  B10  cannot  possibly 
terminate  at  ( 1 )  without  intersecting  B7  and/or  Dg. 

From  II  of  Figure  13,  an  interesting  situation  is  reached  by  letting  (S)  and  (6)  converge,  by 
sequences  of  points,  to  locations  of  (3)  and  (2),  respectively,  such  that  each  polygon  of  the 
associated  sequence  is  SI.  The  limit  element,  however,  shown  in  Figure  15,  is  in  {S}.  There  is 
nothing  contradictory  about  this  with  respect  to  our  previous  remarks.  Figures  (16)  and  (17) 
show  that  S  of  Figure  1 5  can  be  obtained  as  the  limit  of  a  sequence  of  SI  elements,  or  as  the  limit 
of  a  sequence  (Sn(9)}.  Hence,  it  is  still  correct  to  say  the  element  of  Figure  15  is  in  (S),  since  T 
for  it  car.  be  constructed  as  a  simple  polygon  as  shown  in  Figure  17.  And,  it  is  also  correct  to 
maintain  that  if  T  is  SI,  then  there  exists  a  sequence  of  SI  elements  which  converge  to  an  SI  limit 
element  since  by  the  way  T  is  made  there  cannot  exist  a  sequence  (Sn(N)}  converging  to  that  same 
limit  element.  Note  that  S  of  Figure  15  has  a  PAR  at  (4)  which  is  the  reason  for  the  phenomenon 
just  described.  More  will  be  said  on  this  near  the  end  of  this  section. 

We  show  two  more  interesting  examples  in  Figures  18  and  1 9.  The  element  of  Figure  18  is 
in  (S);  i.e.,  it  is  not  SI.  whereas  the  polygon  in  Figure  19.  by  our  characterization,  is  Si. 

SI  polygons  should  not  appear  often  in  practice.  But.  if  the  generation  of  a  polygon  is  not 
under  control  of  the  user,  say  the  nodes  are  computer  assigned,  then  SI  polygons  can  occur.  For 
example,  consider  Figure  20.  It  is  clearly  SI  at  MN(3)  with  self-intersection  point  (3,6).  Note  that 
a  renumbering  1  -*  1.  2  -*  2.  3  -*3,  7  4.  8  -*■  S,  6-»  6, 4->  7.  5  8, 9  -*  9  gives  the  same  regions, 

hut  now  the  renumbered  element  is  in  (S). 


Figure  13.  An  SI  Polygon,  11 
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Figure  20.  An  SI  Poiygoa 


The  basis  for  (P~B)  is  given  by  one  equation,  which  is  a  main  result  of  this  section.  Namely, 
for  any  element  n  of  {11} 

N 

(32)  P<n>  =  W  -  £  P(fli), 

1 

where  P(a,)  is  the  value  of  P  for  the  i,h  exterior  angular  region  of  IT,  and  W  is  a  new  quantity,  which 
we  define  below  and  call  the  winding  number  of  II.  Thus,  there  are  two  basic  steps  here.  The  first 
is  to  evaluate  P(a,)  for  each  i  =  1,  2,  ....  N,  and  the  second  is  to  compute  the  winding  number  of  II. 
The  first  has  been  discussed  extensively  in  the  earlier  sections  and  offers  no  difficulties,  except  that 
there  remains  to  discuss  the  evaluation  of  P(a)  when  a  is  a  PAR.  We  shall  show  below  that  W  can 
be  obtained  by  simply  adding  up  the  angular  measures  A0,  in  radians  of  the  a,  and  dividing  the 
sum  by  2ir.  Thus,  for  an  element  in  {S}  or  {S}  that  is  PO(NO)  W  =  I(-l ),  which  gives  agreement 
of  (32)  with  (24)  ((26))  for  such  elements. 

We  now  need  the  following  definitions  and  additional  notation,  where  n  is  numbered  under 
the  o-opuon: 

A  circuit  of  II  is  a  closed  path  of  11  with  no  self-intersections.  Tlius  a  circuit  is  in  {S>.  and  its 
first  and  last  points  are  located  at  the  same  MN.  Note  that  if  H  is  in  { S>.  then  MN(1 )  is  the  only 
MN  in  the  sense  that  H  is  closed  and  (N  +  1 )  is  at  MN(  1 ). 

A  primary  circuit  of  H»  C,.(U).  is  the  first  circuit  detected  in  tracing  II,  starting  at  (I),  which 
terminates  at  a  self-intersection  of  II, 12  say  between  paths  fk.  k  +  1]  and  |j  *  m,  j  +  m+  li  at 

MN( j ),  where  j  <  k  <  j  +  m,  (see  Figure  1 1  )>  with  CP  «*  (k.  k  +  1 . j  +  m  -  1 ,  j  +  m).  It  contains 

all  other  nodes  (k  +  i)  at  MN(j)  such  that  k  <  k  +  i  <  j  ♦  m.  As  an  example.  C‘P  of  Figure  20  is 
(3.  4,  5.  6);  note  that  k  °  j  °  3,  j  +  in  s  6,  there  are  no  other  nodes  (3  +  i)  at  MN(3)  with 
3  <  3  +  i  <  6.  In  Figure  50.  C‘P  «  (3,  4,  5,  6.  7. 8,  9);  we  have  k  °  3  “  j.  j  +  m  a  9.  all  other  nodes 
(3  +  i)  at  MNP)  with  3  <  3  +  i  <  9  are  included  in  Cp.  Namely,  node  (6)  with  i  ~  3.  i  e,.  CP(II> 
(6,  7.  8.  9).  If  fl  is  an  S  element,  then  CP(n)  £  II. 

The  winding  number  of  a  circuit  C  is  given  by  A0j/2*.  where  At),  is  the  angulav  measure,  m 
radians,  of  the  j,h  exterior  angular  region  of  C.  The  integer  r  denotes  the  number  of  points  speci¬ 
fying  C.  Hie  winding  number  is  + 1  (-1 )  if  C  is  PO(NQ). 

In  order  to  cstablisli  (32),  let  II,  =  M  and  decompose  II  as  follows: 

(c»)  Obtain  CP(H,  V.  Sctial.  Go  to (6). 

(p'i  Find  CP(II,).  Go  to  (5). 

(y)  II  has  been  decomposed  into  a  set  of  disjoint  elements  in  (S3,  (See  Footnote  , 
page  14).  The  decomposition  is  complete.  II  i  s  K.  we  say  ll  lias  or  decomposes 
into  K  primary  circuits. 

(5)  If  (>(»,) 45  II,  go  to  (y).  Otherwise,  delete  CP(II,)  from  II,  (except  for  its  last  node), 
call  the  result  II,,, .  set  i  ♦  I  =  i.  and  go  to  (jJ). 

,JCare  roust  be  taken  when  PAR(s)  ate  involved.  More  later.  See  footnote  1 3. 
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For  example,  for  the  SI  polygon  of  Figure  20,  we  have 

rij  =  (i,  2, ...,  9,  io),  Cpdi,)*'?.  4, 5,  6),  n2  =  n,  -cP(n1)  =  (i,  2,6, 7,8,9, 10) 
Cp(n2)  =  n2.  Thus  II  =  Cpfllj )  U  Cp(II2),  where  U  denotes  union. 

The  decomposition  of  II  into  primary  circuits  can  always  be  carried  out.  Indeed,  by  a  slight 
modification  of  a  proof  by  Knopp  [4,  page  15] ,  one  can  state  that  any  polygon  can  be  decomposed 
into  a  finite  set  of  disjoint  elements  in  {S>.  Knopp’s  proof  is  constructive,  and,  if  it  is  followed,  as 
described  in  Appendix  A  for  (P  -  A),  the  polygon  II  is  decomposed  into  a  set  {§}  =  (S1,  S2,  S3, ..., 

SJ)  of  simple  polygons  and  a  set  { L}  =  (L1,  L2 . L°)  wheie  each  L1  denotes  an  overlapping  line 

segment  (a  PAR).  Now  CP(1T t )  is  made  up  of  the  union  of  a  subset  of  {§}  and  a  subset  of  {L}. 
Deleting  these  elements  from  {S}  U  {T}  leaves  II2.  Then  CP(II2)  is  found,  for  II2,  in  the  same  way 
from  {§}  U  {L>  -CpOI,).  For  example  in  Figure  22  we  have  L1  =  (3,  4,  5),  S1  =  (1,  2,  3,  6)  and 
cP(n1)  =  L1  us1  =  n,  with  Wj  =  w  =  i. 


In  general, 


K 

(33)  n  =  (J  CpOI;),  1  <  K  <  N, 

t 

where  the  Cpfl!;)  are  disjoint,  (See  Footnote  9,  page  14).  Hence 

K 

(34)  P(II)  =  T  P[CP(ni)]  , 

l 

where,  using  (24)  and  (26), 

f  1  if  CpOlj)  is  PO 

P(*in).  W,  = 

{  -1  if  CpOI,)  is  NO, 

with  <2jj  denoting  the  jth  of  Kj  exterior  angular  regions  of  Cp(flt).  Then,  substituting  (35)  into  (34) 
gives 

K  K 

(36)  P(H)  Wi  ~L  E  P<fl»n)* 

I’l  iBl  0»l 


K| 

(35)  P[Cp(n,)]  =  Wj  -  2 


The  winding  number  of  II,  W,  is  defined  to  be 

K 

(37)  Ws^W,. 

i*i 
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Now  let 


(38) 


N 

£2  =  A 0j,  -ir  <  A0;  <  7r,13 

i=l 


where  AO:  has  the  usual  meaning,  i.c..  it  denotes  the  angular  measure,  in  radians,  of  the  exterior 
angular  region  of  n.  To  establish  (32)  and  that  W  can  be  expressed  in  terms  of  £2,  we  need  to 
show  that 


(39) 


r  S2/2-?  a  w 


< 


N  K  K, 


E  p<3»> ■  E  E 

k  =  l  1  =  1  n=l 


We  proceed  to  present  the  elements  of  a  proof  for  (39).  The  argument  goes  as  follows: 
Suppose  Cp(ilj )  =  (j,  j  +  1,  ....  j  +  m),  so  that  n2  =  ( 1,  2,  ...,  j  -  1,  j,  j  +  m  +  1, ...»  N  +  1),  where 
n2  =  n,  -  CP(II|).  Denote  the  exterior  angles  of  CP(I1])  and  n2  at  (j)  by  a,  j  and  a2  j, 
respectively.  Denote  their  corresponding  angular  measures  by  f  and  >.  Also,  let  A^  and  A 0V 
t  ~  j  +  m,  denote  the  measures  of  a}  and  a,  of  n. 

With  the  aid  of  Figure  2 1  below,  we  have  for  one  particular  situation, 

(40)  J  +  X  =  A0j  +  Afl, . 

Then  from  geometrical  arguments,  we  must  have  also 

(41)  flj.,  U  a2 ,j  =  cfj  u  aj,m.  and  P(aKl)  +  P(«2.))  =  P(ap  + 


,JTo  remain  consistent  with  our  characterization  of  S  and  SI  elements  front  the  T -construction,  wc  need  IAfl(|  <  n, 
but  computationally  wc  always  have  -s  <  A0j  <  tr  as  a  result  of  using  the  4*quadranl  arctangent  routine  to 
compute  AOj. 
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In  fact  (40)  and  (41)  hold  for  any  other  configuration  with  a  self-intersection  at  (j,  j  +  m).  There  are 
five  other  basically  different  configurations  to  check.  This  has  been  done;  the  details  are  not  given. 

Assume  now  that  II  has  only  one  self-intersection.  Then  by  the  relations  given  in  (40) 
and  (41),  we  see  that  the  only  angular  regions  affected  are  at  (j)  and  (j  +  m),  where  the  intersection 
takes  place.  Therefore  with  (34)  and  (35) 

(42)  P(n)  =  PiCpOij)]  +  P(n2),  (w  =  Wj  +w2), 


Wi  - 

j+m  - 1 

p(ou>+  E  ffek) 

+ 

* 

1 

'j-1 

y  P(ak)  +  P(a2fj)  + 

N 

E  pk> 

k=j+  1 

_k=  1 

k=j+m  +  l 

Using  (41)  in  (42)  we  have 


(43) 


P(II)  =  W 


E 


k  =  l 


P fefc)  • 


It  remains  to  show  the  first  equation  of  (39)  in  the  case  of  one  self-intersection.  We  have  the 
winding  numbers  Wt  and  W2,  for  Cpdlj )  and  Il2)  respectively,  given  by 


(44) 


2ttWj 


j*m-t  j-t  N 

f  +  £  A0k ,  2irW2  =  X  +  Y  A0k  +  T  A0k, 
j-H  1  j+m+l 


where  II  with  only  one  intersection,  implies  that  Il2  is  a  circuit.  Consequently,  using  (37)  and  (40) 


(45) 


W  =  W,  +  W2 


JL 

2jt 


?  +  X  +  £ 

kBl 

k*J,J+m 


12/  2*r. 


To  treat  the  case  where  II  decomposes  into  K(>2)  primary  circuits,  an  induction  argument 
can  be  used.  Assume  (32)  holds  for  all  elements  n  with  no  more  than  K  -  1  primary  circuits.  The 
essence  of  a  proof  that  (32)  holds  for  polygons  with  (decomposable  into)  K  primary  circuits  is 
obtained  from  the  argument  above  for  K  "  2. 


Let  n  have  K  primary  circuits.  Then  by  the  decomposition  procedure  described  above, 

n  =  Cp(n,)  u  n3. 

where  II2  can  be  decomposed  into  K  -  l  primary  circuits  with  winding  numbers  W3 ,  WJ(  ...,WK. 
But  (32),  by  the  induction  hypothesis,  holds  for  II2.  Therefore,  the  remainder  of  the  argument 
goes  as  above  for  K  °  2,  where  W2  is  replaced  by  Wt  in  (42),  (43),  (44)  and  (45). 


On  pages  10,  12,  we  remarked  that  the  evaluation  of  P(crk)  requires  seme  precaution,  when 
ak  is  a  PAR,  ie.,  I  A0k  1  =  i.  The  basic  difficulty  is  tliat  the  sign  of  A0k  cannot  be  determined  from 
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the  arctangent  subroutine,  because  it  always  gives  it  for  a  PAR.  Thus  to  determine  the  sign  requires 
some  additional  analysis.  This  analysis  amounts  to  correctly  choosing  the  direction  Tk+I  should  take 
in  the  construction  of  T  (see  page  1 5).  For  example,  an  attempt  to  construct  T  for  S  of  Figure  1 5 
by  choosing  T5,  as  shown  in  Figure  16,  would  not  be  right.  This  amounts  to  choosing  the  sign  of 
A04  incorrectly  since  elements  typical  of  the  one  in  Figure  16  are  SI.  The  correct  direction  for 
T5  is  shown  in  Figure  17  which  indicates  A 0A  =  -it  rather  than  7r. 

We  can  further  elucidate  the  difficulty  with  the  aid  of  Figures  22  and  23.  The  ele¬ 
ments  in  both  figures  are  in  S.  Call  them  S2  and  S3,  respectively.  Both  elements  have  a  PAR 
at  (4).  In  Figure  22,  the  value  for  A04  of  S2  is  tt  and  for  S3,  A04  =  -it.  We  see  for  S2  that 
1  -  P(S2)  =  P(tfj)  +  P(a2)  +  P(a),  and  P(a)  -  P(a4)  +  P(a3)  +  P(tfj).  Thus  any  routine,  such  as 
VALR-7,  designed  to  compute  P  for  elements  in  {S},  but  not  for  SI  polygons,  would  give  the 
correct  result  for  S2 .  For  S3  however,  the  situation  is  different.  The  arctangent  subroutine  would 
give  A04  =  it  which  is  incorrect  for  S3  to  be  in  {S}  and  consequently  P(a3 )  +  P(a4 )  +  P(a5 )  =£  P(a) 
and  the  result  from  VALR-7  for  P(S3)  would  be  wrong. 

The  subroutine  VALR-7  forms  a  part  of  our  preferred  program  package,  because  of  its  slightly 
superior  speed  over  VALR-2  in  computing  P  for  elements  in  {S}.  So,  in  order  to  also  use  VALR-7 
for  elements  in  fS>,  rather  than  include  the  additional  steps  in  the  program  to  detennine  AO  cor¬ 
rectly  for  PAR(s),  a  routine  SORT  III  was  designed  which  pre-processes  S  by  deleting  from  its 
specifying  array  V  all  SAR(s)  (see  page  13).  The  result  is  an  element  in  {S>,  say  S,  such  that 
P(S)  =  P(S),  since  SAR(s)  do  not  affect  the  value  of  P(S)  by  their  removal. 

In  the  case  of  VALR-2.  which  is  based  on  (P-B),  the  winding  number  rectifies  any  wrong 
choice  for  A 0  at  a  PAR.  For  S2,  in  Figure  22,  A04  is  computed  correctly ;  therefore,  W  =  1  for  S2. 
For  S3  however,  A04  is  computed  incorrectly  as  noted  above,  so  that  from  Figure  23,  we  have 

P(u4)  +  P(a3)  +  P(tfj)  =  1  +  P(fl), 

where  P(u)  is  the  correct  value  to  add  to  P(at )  +  P(a2)  to  get  P(S3).  Now  by  (32) 

P(Sj)  =  W  -  |P{a,)  +  P(a2)  +  1  +  P(a))  . 
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Figure  23.  A  Polygon  S3  in  {S}, 
(Ad 4  incorrect),  W  =  2 


But  note  that  the  computation  of  W  using  (38)  with  the  arctangent  routine  is  two,  since  Cpd^)  = 
(3,  4,  5),  CP(n2)  =  (1,  2,  3,  6)  and  Wl  +  W2  =  1  +  1  =  2.  Hence  P(S3)  is  given  correctly  by 

P(S3)  =  1  -  [P (a,)  +  P(a2)  +  P(a)l . 

Hence,  an  element  with  PAR(s)  may  be  classified  by  the  T-construction  as  in  (S>  and  yet  its 
winding  number  from  (38),  using  the  arctangent  routine,  (See  Footnote  13)  will  not  necessarily  be 
i  1 .  Thus  computationally  it  must  be  treated  as  SI. 

It  is  to  be  noted  that  the  value  of  P  for  a  PAR  requires  an  erfc  function  computation.  For 
example,  in  Figure  30,  W  =  6  and  8  erfc  functions  are  needed.  Hence,  it  may  be  worth  using 
SORT  III  also  with  VALR-2  to  eliminate  SAR(s),  (See  flow  chart  for  P-2,  page  40). 

We  conclude  this  section  by  giving  two  more  examples  of  decomposing  II  into  primary  circuits, 
and  some  remarks  on  numbering  II  other  than  by  the  ooption. 

Consider  the  polygon  of  Figure  24.  We  shall  show  that  W  »  3.  Note,  the  first  self-intersection 
occurs  at  (4,7): 

Cp(ll,)  =  (4,5,6, 7),  W,  =  1,  Ilj  -  (1,2,  3,  7,  8,  9,  10,  II,  12,  13,  14,  15). 
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The  first  SI  point  of  II2  is  at  (3, 1 1);  by  the  preceding  discussion, 

CP(n2)  =  (3,  7,  8,9,  10,  11),  W2  =  1,  n3  =  (1,2, 11,12, 13, 14,  15), 
where  CP(I12)  is  PO  and  in  {S}.  With  no  remaining  self-intersections,  we  have 

cP(ii3)  =  n3,  w3  =  i  . 


Hence  W  =  3. 


For  Figure  25  (note  the  a-numbering)  the  decomposition  of  II  is  as  follows: 


CpOl,)  =  (6,  7,  8,  9,  10,  11), 

Wi  =  -1 , 

CP(n2)  =  (4,5,  11,  12), 

w2  =  1, 

Cp(II3)  =  (2,3,  12,13,14,15), 

W3=-l, 

Cp(II4)  =  n4,  W4  =  1 . 

n2  =  (1,2,3,4,5,11,12,13,14,15,16,17) 
n3  =  (1,2,3,  12,  13,  14,  15,  16,  17) 

II4  =  (1,  15,  16,  17) 


Hence  W  =  0. 

In  the  actual  computation  of  P(II),  it  is  often  not  necessary  to  specify  n  by  an  a-numbering 
(see  page  14).  However,  one  should  not  assume,  for  example,  that  all  points  for  which  AO  55  0  can  be 
dropped  when  using  (P~A),  although  this  is  permissible  in  (P-B).  In  Figure  55,  if  that  element 
were  treated  as  SI,  an  additional  node  at  MN ( 1 )  would  be  required,  following  (6),  for  (P-A),  using 
SORT  l,  to  give  the  correct  result.  This  will  be  clarified  in  Appendix  A  (See  page  A-7). 

In  Figure  25,  six  points  can  be  dropped  for  computation  purposes  using  (P-B),  namely  (2). 
(4),  (6),  (11),  (12),  (15).  The  reduced  numbering  scheme  as  shown  in  Figure  26  is  called  a 
^-numbering  scheme,  cr  simply  a  Option. 

An  additional  short  discussion  on  the  superiority  of  (P-B)  over  (P-A)  is  given  on  page  A- 1 5. 


Figure  25.  SI  Polygon  with 
a-Option.  W  =  0. 
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Figure  26.  SI  Polygon  with 
0-Option.  W  =  0. 


V.  DISCUSSION  OF  COMPUTER  PROGRAM  B  (FLOW  CHARTS  INCLUDED) 

The  program  package  described  here,  call  it  Program  B,  contains  five  subprograms,  each  in  sub¬ 
routine  format,  P-2,  VALR-2,  SORT  III,  VALR-7,  SMP-7.  The  second,  fourth,  and  fifth  can  be  used 
independently;  the  first  serves  as  a  master  routine.  The  last  three  are  also  used  in  the  program  of 
Appendix  A,  called  Program  A.  We  shall  discuss  each  subprogram,  and  point  out  how  each  is  used 
in  connection  with  the  others.  VALR-7  has  much  in  common  with  VALR-2;  a  detailed  discussion 
on  it  is  not  given.  However,  the  essential  differences  between  it  and  VALR-2  are  noted.  Keep  in 
mind  that  VALR-7  is  on  the  average  the  slightly  more  efficient  of  the  two  for  computing  P  for  ele¬ 
ments  in  {S},  but  VALR-2  alone,  which  is  based  on  (P-B),  is  capable  of  determining  P  for  any 
element  of  {II}.14  Flow  charts  for  the  first  four  subprograms  are  included  at  the  end  of  this  section, 
pages  40-45.  They  will  be  used  to  aid  in  the  discussion.  No  flow  chart  is  given  for  SMP-7,  which  is 
used  to  compute  A. 

The  given  polygon,  call  it  n,  for  which  P  is  wanted,  is  specified  by  its  nodes  at  points 
(xk,  yk),  k  =  1,2, ....  N.  The  call  line  of  each  routine  identifies  these  quantities  by  x,  y,  N. 

The  parameter  IOP  appears  in  the  call  line  of  P-2,  VALR-2,  VALR-7.  It  specifies  the  approxi¬ 
mate  accuracy  to  which  P(H)  is  computed.  The  user  assigns  IOP  =  1, 2,  or  3,  so  that  P  for  each 
angular  region  of  II  is  computed  with  3, 6,  or  9-decimal-digit  accuracy,  respectively. 

The  parameter  ICV  appears  in  the  call  line  of  P-2.  With  this  parameter,  the  user  can  specify,  for 
maximum  efficiency  of  computation,  various  combinations  of  the  above  routines  or  subprograms. 
The  flow  chart  for  P-2,  page  40,  summarizes  the  combinations  one  may  choose.  If  II  is  in  {S},  then 
the  user  should  set  ICV  =  0  and  P-2  would  call  VALR-7  to  compute  Pdl).  When  N  =  1  is  specified, 
the  normal  probability,  regardless  of  the  ICV  value,  is  given  foi  the  angular  region  a  which  is 
specified  by  three  points.  If  the  user  is  uncertain  about  the  class  to  which  H  belongs,  ICV  >  0 
should  be  used.  Then  P-2  calls  VALR-2  to  find  P(ll)  where  fl  can  be  in  {II}.  If  II  has  tr-angulur 
regions,  PAR(s)  (see  page  12)  and  II  is  not  SI,  then  it  may  be  more  efficient  to  use  SORT  III  with 
VALR-7  rather  than  VALR-2  alone.  This  can  be  done  by  setting  ICV  =  -2.  If  II  has  self-in  terscctions 
as  well  as  PAR(s),  then  VALR-2  with  SORT  111  may  be  more  efficient  ‘han  VALR-2  alone  since 

MThc  package  above  gives  an  improvement  in  efficiency  over  using  VALR-2  alone. 
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VALR-2  makes  an  erfc  function  computation  for  each  PAR  of  II;  SORT  III  deletes  such  regions 
before  VALR-2  is  called.  This  combination  can  be  called  by  setting  ICY  <  0,  but  not  equal  to  -2. 
In  using  SORT  III,  N  may  be  reduced  below  3,  in  which  case  P  =  0,  A  =  0  as  shown  in  boxes  9 
and  13  of  the  flow  chart  for  P-2,  page  40. 

Later,  in  the  discussion  of  VALR-2,  the  reason  why  VALR-7  can  be  slightly  more  efficient 
will  be  given.  But,  VALR-7  can  give  grossly  wrong  values  of  P  if  it  is  used  for  SI  polygons  or 
without  SOR  T  III  for  an  S  element  with  SAR( s),  (See  page  31). 

A  by-product  of  VALR-2  and  a  necessary  quantity  for  VALR-7  is  A,  which  is  given  by 
1  N 

(46)  A  =  j  £  xi(yi+i yo  s  Vn>  yw+i  s  Yi  *  N  >  3. 

l 

For  VALR-7,  A  is  computed  by  the  subroutine  SMP-7.  It  is  used  in  VALR-7  to  determine  the 
orientation  of  IT,  when  IT  is  in  {S}  or  (S).  The  sign  of  A  determines  whether  P(I1)  in  VALR-7  is 
evaluated  by  (24)  or  (26),  (See  page  9). 

No  flow  chart  is  given  for  SMP-7,  but  a  listing  of  the  program  is  given  in  Appendix  F, 
page  F-37.  A  derivation  of  (46)  is  given  in  Appendix  D,  where  it  is  shown  that  |  A 1  when  properly 
interpreted  gives  the  area  of  II. 

For  VALR-2,  A  is  computed  within  VaLR-2  itself.  In  VALR-7  it  plays  a  crucial  role  as 
evidenced  by  (24)  and  (26).  The  winding  number  of  IT,  W,  plays  a  similar  but  more  complicated 
role  for  VALR-2  as  (32)  shows.  In  the  previous  section,  it  was  shown  that  W  is  computed  from 
2,N  A0,/2ir.  where  A0(  is  the  angular  measure  of  a s,  the  exterior  angular  region  of  II  at  (i).  The 
winding  number  in  addition  to  Pdl),  A(I1),  and  IND  make  up  the  output  of  VALR-2. 

The  error  indicator  IND  is  used  in  both  VALR-2  and  VALR-7.  Its  normal  setting  is  zero.  If 
VALR-2  is  used  to  find  P  for  an  angular  region  (N  =  1),  and  the  x,  y  input  contains  a  SDP,  then  IND 
is  set  to  one  and  P  is  set  to  the  absurd  result  of  5.  If  IND  is  set  to  two  in  VALR-7,  or  VALR-2.  it 
means  a  PAR  was  encountered  in  evaluating  P(H),  and  the  result  for  P(I1)  from  VALR-7  is  not  to 
be  trusted,  whereas  the  corresponding  result  from  VALR-2  is  good.  If  IND  =  3,  then  a  direct  exit  is 
taken  since  N  has  been  set  as  smaller  than  one  or  equal  to  two.  In  either  case  such  an  assignment 
is  not  acceptable  to  VALR-2  or  VALR-7.  See  boxes  7  and  l  in  flow  charts  of  VALR-2  and  VALR-7, 
respectively. 

For  easy  reference,  the  above  programs  are  numbered  accordingly:  P-2  1,  VALR-2  «+  2, 

SORT  III  <*  3,  VALR-7  4,  SMP-7  <♦  8. 

We  proceed  with  a  more  detailed  discussion  of  2.  We  refer  to  the  n,h  box  of  the  flow  cltart 
for  2  by  2-n. 
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Although  4  should  be  used  when  N  =  1,  this  case  is  also  handled  by  2.  When  N  =  1,  P(a) 
as  computed  by  2  (or  4)  always  gives  the  normal  probability  for  a ;  a  negative  value  is  never 
found  for  P (a)  in  this  case.  Three  points  are  necessary  to  define  a,  with  (1)  always  referring  to  the 
vertex  of  a  with  points  (2)  and  (3)  given  in  counterclockwise  order  with  0  <  A0  <  2n.  Notice  that 
this  differs  from  A0  for  an  exterior  angular  region  of  a  polygon,  where  |  A0|  <  ir.  In  Figures  27  and 
28,  the  assignments  of  3(x,  y)  coordinates  are  shown  for  two  different  angular  regions  when  N  =  1. 

The  sensing  for  N  =  1  is  carried  out  at  2-7  (flow  chart  for  VALR-2,  box  7). 15  If  *//,  defined  in 
2-15  (see  also  pages  28  and  39),  is  nonnegative,  we  have  an  element  of  type  shown  in  Figure  27, 
0  <  Ad  <  n.  If  \p  <  0,  Figure  28  represents  a  typical  case,  tt  <  A6  <  2ir.  If  a  0,  then  a  can  have 
a  vertex  angle  near  0,  ir,  or  2tt  radians.  Here  the  user  must  be  careful,  because  if  rounding  error 
should  .interchange  (2)  and  (3),  a  wrong  result  can  be  given  for  P(a),  (See  page  8).  The  error 
indicator  IND  is  not  used  for  this  situation.  The  other  boxes  used  only  for  N  =  1  are  2-57,  and  2-79. 

If  N  >  3,  2  treats  any  polygon  II  by  finding  P(ak),  k  =  1,  2,  ...,  N,  for  each  exterior  angular 
region  ak  of  IT.  The  analysis  for  computing  P(ak),  using  (13)  and  (17),  was  given  in  Sections  II  and 
III.  The  A0k  term  for  ak  (see  Figures  6,  9,  pages  5,  10)  is  computed  at  2-36  using  a  4-quadrant 
arctangent  routine  which  gives  values  in  the  half-closed  interval  ( — 7r,  7r] .  The  sum  of  the  A0k  is 
accumulated  in  $2  at  2-27.  Individual  terms  for  A,  as  given  in  (46),  are  computed  at  2-22  (k=  1) 
and  2-68  and  accumulated  in  A.  The  sum  in  (13)  is  computed  at  2-35  and  the  loop  2-34,  54.  The 
am_j  in  2-34  are  the  Chebyshev  coefficients  for  erfc(u)/z(u)  as  they  appear  in  (13).  They  are 
tabulated  in  Appendix  E  for  the  three  IOP  settings,  with  an  additional  set  for  1 2-decimal  digits  of 
accuracy  which  is  not  included  in  the  program.  The  value  of  P(ak)  is  given  by  I  in  2-47,  where  L 
refers  to  the  crfc  function  contributions  from  using  (17).  Note  if  gj  and  g2  from  2-19  are  non¬ 
negative  2-29, 30,  then  (1 7)  is  not  needed  and  L  °  0,  2-23,  2-37. 


15 The  boxes  in  each  (low  diart  are  numbered,  usually  at  the  upper  right-hand  corner. 


For  efficiency,  yet  retaining  specified  accuracy,  if  R  is  sufficiently  small  or  large,  then  ( 1 3)  is 
not  used.  This  is  reflected  through  the  sensings  at  2-18,  1 2,  28.  If  R2/2  <  a,  then  I  =  Ad/2ir  and 
if  R2/2  <  a2,  then  I  =  (Ad/2ir)  -G,  where  G  is  computed  at  2-10  from  (18).  If  R  >  R,  then  I  =  L, 
2-28,  37.  At  this  point,  VALR-7  can  be  more  efficient  than  2.  In4(VALR-7),  A0k  is  not  computed 
until  needed  (see  4-16,  34)  so  that  when  (13)  is  not  used  an  arctangent  computation  is  saved.  In 
2,  however,  A 0k  is  computed  regardless  of  whether  (13)  is  used  or  not,  because  it  is  needed  to  find 
the  winding  number  W  =  Sl/lit.  Thus,  for  a  polygon  with  n  vertices  of  S  located  such  that  R  >  R 
for  each  of  them,  2-28,  means  n  more  arctangents  would  be  computed  by  2  than  by  4.  The  param¬ 
eters  cq,  of2,  R2/2  as  well  as  a3  and  a4,  which  are  discussed  below,  are  given  in  Appendix  E.  They 
depend  on  the  setting  of  IOP. 

In  2-19,  the  rotation  of  axes  for  ak  is  done,  as  discussed  on  page  4. 

If  s,  given  by 


(47)  s  s  2^/D,D2  =  sin  A016 

is  sufficiently  small  in  absolute  value,  2-20,  then  ak  subtends  an  angle  near  0  or  ±rr.  When  it  is 
sufficiently  close  to  zero,  with  <j>>  0,  2-24,  and  |  A0|  <  7 (-14),  2-14,  then  0  — ►  I,  2-11.  If  there  is  a 
no  at  2-14  then  I  ss  0  only  if  >0,2-13.  If  ga  <  0,  with  A0  ~  0,  it  may  happen  that  ak  contains 
the  origin,  in  which  case  P(ak)  is  not  near  zero  for  sufficiently  large  R. 

If  Isl  <cv3/4,  2-20,  with  |A0|  -  tt(0  < 0,  2-24),  and  if  |s|  «7(-14),  IND  is  set  to  two,  2-21. 
There  is  nothing  to  be  concerned  about  here.  IND  is  simply  used  to  alert  the  user  that  a  PAR  has 
been  encountered.  Recall  that  in  4  if  this  occurs,  P(II)  is  not  to  be  trusted,  (This  is  discussed 
further  on  page  A-3  of  Appendix  A).  Now  if  \J/  <  0,  2-6,  then  P(ak )  is  given  by  1  of  2-4,  and  if 
i//  >  0,  then  P(ak)  is  given  by  l  of  2-5.  Boxes  4  and  5  are  where  the  erfe  function  computations  are 
made  in  the  program  for  a  PAR.  Note  the  choice  is  made  such  that  if  ^  =  0  it  is  assumed  A 0k  =  ir. 
Keep  in  mind  that  if  111  I  denotes  the  norma!  distance  from  an  extended  straight  line,  intersecting 
the  nonnegative  x-axis,  to  the  origin,  then  1/2  erfc  (h/\/2)  with  h  >  0«0)  gives  the  probability  over 
the  half-plane  to  the  right  (left)  of  the  line. 

Assume  now  that  the  program  moves  from  2-20  or  2-13  to  2-12  and  from  there  to  2-29.  Then, 
in  the  next  set  of  boxes,  starting  at  2-29,  the  necessary  erfc  computations,  required  by  (17),  or 
approximations  to  them  are  made  and  stored  in  L.  As  mentioned  above  if  g,.  g2  >0,  then  (17)  is 
not  used,  L  is  set  to  zero,  and  control  is  transferred  to  2-28. 


We  use  the  notation 


(48) 


H(h)  =  erf(h)  =  1  -  erfc(h)  -  1  -  E(h),  (see  (8)), 
E(-h)  -  1  +  0(h)  =  2  -  E(h) 


'^Variables  appearing  in  the  flow  charts  arc  defined,  or  cross-referenced  to  the  text,  on  page  39,  which  precedes  die 
flow  diarts  of  this  section. 
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If  gj  <  0,  g2  >0,  then  -gj  ->  gj,  -hj  -*■  hj,  2-29,  38  and  2-32. 

If  ^<0,  then  A6  -ir^  Ad,  2-25,  and  L  =  -j  E(ht),  2-45. 


If  \jj  >  0,  then  A d  +  it  ^  A 6,  2-41,  and  L  =  -j  E(-ht). 

If  gj  >  0,  g2  <  0,  then  -g2  -*■  g2,  -h2  -*•  h2,  2-29,  30,  39,  and  ~it  -+  \p,  2-32. 

If  *  <  0,  then  A6  ~n  -*  Ad,  2-25,  and  L  =  j  E(-h2),  2-31. 

If  \p  >  0,  then  Ad  +  ir  -*  AS,  2-41  and  L  =  -■j  E(h2) 

If  gj  <  0,  g2  <  0,  then  -g1  ->  gj,  -hj  -*•  hlt  2-29,  38,  and  -g2  -*  g2,  -h2  -*  h2,  243,  and 
L  =-  (E(hj)  -  E(h2)],  2-52. 

An  approximation  for  E(h)  or  E(h)  is  used  for  efficiency  of  computations,  if  1  h  |  is  sufficiently 
small,  i.e.,  |h|  <  a4.  The  sensings  on  this  inequality  occur  in  2-39, 43, 44, 49,  50.  In  each  case,  if 
the  inequality  is  satisfied,  E(h)  is  replaced  by  (2/-n/tt )h.  The  parameter  a4  depends  on  10P,  and  is 
determined,  to  tetain  the  accuracy  specified,  by  using  the  mean  value  theorem  on  E(h).  Indeed, 

(-t2)  dt  =  — 

\ /» 

Retaining  the  first  term,  the  error  e(h)  is  bounded  accordingly: 

(50)  |o(h)l  =-^|4$2-2|  ^  exp(-f2)  <%. 

Thus 


h^ 


h  +  (4*2-2)  ~  exp  (~£2) 


{e(0,  h). 


(49)  E(h)  a  J  exp 


(51) 

with  (49)  and  (50),  implies 

(52) 


ih|  <  a4  s  e 


1/3 


E(h)  -  ~p  h 
V* 


<  e/2 . 


where  e  and  oi4  are  given  in  Appendix  E  as  a  function  of  the  10P  setting.  For  example,  if  (51)  is 
satisfied  for  h,  and  h2.  with  ,  g2  <  0,  then  L  is  computed  from  248  rather  than  2-52  with  an 
error  no  larger  than  e. 


After  L  has  been  determined,  control  proceeds  to  2-28  to  check  if  R  is  sufficiently  large  so 
that  the  calculation  of  (13)  can  be  omitted,  or  if  the  second  term  on  the  right-hand  side  of  (17) 
can  be  dropped. 

Consider  the  angular  region  a  shown  in  Figure  29,  which  corresponds  to  case  JIT]  in  Figure  5 
of  [2,  page  14].  The  quantity  P(a)  is  found  using  (17),  i.e.,  P(a)  =  l/2[E(h2)  -  E (h: )]  +  P(a), 
where,  at  2-19,  gj,  g2  <  0,  hj  >  G,  h2  <  0.  Now  if  R  is  sufficiently  large,  say  >R,  2-28,  thenP(a ) 
is  negligible  and  its  computation  by  (12),  2-35,  34,  47,  can  be  bypassed  by  proceeding  directly 
to  2-37. 

In  case  R  <  R,  2-28,  the  quantity  (13)  is  computed  from  2-35,  34,  47.  Then  P(ak)  is  obtained, 
as  shown  in  247,  with  the  result  stored  in  I.  Recall  that,  for  efficiency,  P (ak)  may  be  computed  in 
various  other  ways  as  indicated  at  2-2, 4,  5,  11,  37. 

Control  is  now  transferred  to  2-55  to  determine  if  N  angular  regions  have  been  processed.  As 
each  P(ak)  is  computed,  it  is  subtracted  from  P,  which  is  initially  set  to  zero,  2-66.  If  the  answer  is 
no  at  2-55,  the  quantities  w,  z,  Dt  at  2-59  and  u,  v,  D2  are  updated  at  2-64  to  be  used  for  the  next 
angular  region.  Then  P  —  I  -*■  P  is  carried  out  at  2-66,  as  noted  above,  and  x(yk+1  -  y)  +  A  -*■  A, 
2-68  also  noted  earlier.  Control  is  then  returned  to  2-26. 

When  [-2J1  P(ak)]  has  been  computed,  i.e.,  k  =  N,  2-55,  56, 60,  then  W  =  S2/2jt  is  computed, 
2-61.  The  quantity  P(fl)  is  then  found  from 

N 

(53)  P(n>  =  W  -  £  P(akt 

i 

at  2-63  or  2-67  (see  (32)),  with  W  now  an  integer  variable. 

The  remaining  boxes  of  the  flow  chart  2  to  discuss  deal  with  the  handling  of  $DP( successive 
duplicate  points).  Whenever  two  successive  nodes  of  II  occur  at  the  same  xy-point,  one  of  them  is 
ignored  in  computing  P(U).  This  feature  is  not  contained  in  VALR-7,  since  it  is  handled  by 
SORT  III.  It  is  included  in  VALR-2  so  that  this  subroutine,  entirely  on  its  own,  can  find  P(I1)  for 
any  H  in  {II}. 


Figure  29.  Shows  Angular  Regions  a  and  a 
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If  Dj  <  o>2,  where  to  =  7(-14),  k  =  1,  2-16,  then  (1)  and  (N)  are  SDP  (within  to)  and  control 
is  transferred  to  2-75.  If  initially  N  =  1 ,  then  IND  is  set  to  1  and  P  to  5  since  the  angular  region  a  is 
not  defined,  2-79.  If  N  >  3,  then  N  =  N  -  1  and  if  N  is  reduced  to  two,  P  is  set  to  zero  for  IT  with 
an  exit  2-77,  80.  If  N  is  not  reduced  to  two,  control  is  returned  to  2-8,  where  new  values  of  w  and  z 
are  computed  and  Dj  is  checked  again,  2-16.  If  (1)  and  (N)  are  not  SDP  then  D2  <  <o2  is  checked 
at  2-17.  If  this  inequality  is  satisfied,  then  (1)  and  (2)  are  SDP  and  control  is  transferred  to  2-81. 
If  N  =  1,  again  the  angular  region  is  not  well  defined,  IND  =  1,  P  =  5  and  EXIT,  2-79.  If  N  >  3, 
then  k  +  1  =  k  and  new  u  and  v  are  computed,  2-78,  72,  and  the  inequality  D2  <  cc2  is  checked, 
2-73.  If  it  is  satisfied  return  to  2-78,  increase  k  by  one,  and  repeat  2-72,  73.  A  no  eventually 
must  occur  at  2-73,  because  at  k  =  N,  points  ( 1 )  and  (N)  are  checked,  but  these  cannot  be  SDP  since 
for  k  =  1  they  were  checked  at  2-16.  With  a  no  at  2-73,  is  k  =  N  -  1?  If  the  answer  is  yes,  then 
points  (2),  (3),  ...,  (N  -  1)  are  each,  with  (1),  SDP,  so  P  =  0,  and  EXIT  is  made,  2-74,  83,  80.  If 
the  answer  is  no  at  2-74,  then  (N)  and  (1),  (1)  and  (k  +  1)  in  the  array  specifying  II,  are  not  SDP; 
control  is  returned  to  2-22,  and  2  proceeds  to  compute  P(a)).  Computation  of  new  u,  v,  D2  quanti¬ 
ties  and  updating  of  w,  z,  Dj  at  2-59,  64  for  k  >  1;  also  includes  a  check  to  see  if  (k)and  (k  +  1) 
are  SDP. 

The  barred  x  and  y,  2-1,  59,  64,  68,  and  the  y,  2-59,  68,  are  used  so  that  once  two  SDP  are 
found  testing  for  more  such  points  at  the  same  k  can  be  done  with  reference  to  the  same  point, 
namely  (x,  y). 

We  next  look  at  SORT  III.  Its  function  is  to  remove  points  from  V,  the  xy-array  which 
specifies  IT,  when  either  of  the  following  conditions  hold  in  considering  ak 

(A)  Either  k  or  k  +  l  (see  page  1 4)  has  a  length  no  larger  than  to  =  7  (- 14). 

(U)  The  angle  AQk  subtended  by  ak.  satisfies  one  of  the  inequalities  ir  -  to  <  A0k  <  tr, 

-v  <  <  to  -  v. 

If  (A)  holds,  we  say  (k  -  I)  and  (k)  or  (k)  and  (k  +  l)  are  successive  duplicate  points  (SDP).  If  (B) 
holds,  we  say  (k  -  I),  (k),  and  (k  +  l)  form  a  PAR.  In  either  case,  we  say  ak  is  a  singular  angular 
region.  (SAR).n  because  ir  (A)  holds  I)2  a  w-  +  t}  or  D2  -  u2  +  v2  is  essentially  zero,  yet  they 

must  be  used  as  divisors  in  2-19  (or  4-24);  if  (B)  holds  then,  because  the  arctangent  subroutine 

gives  values  on  (“*.  tr| ,  4  cannot  determine  whether  AO  s  it  o r  f— nr>  for  a  PAR  (this  is  discussed  on 
page  22).  For  2,  (B)  holds  no  difficulty  as  explained  on  page  22.  If  neither  (A)  nor  (B)  holds,  we 
say  ak  is  well  defined  (WD). 

Given  an  array  V  made  up  of  the  ordered  set  of  N  xy-coordinalcs  which  define  11.  SDP  are 
climi;  ted  by  dropping  one  of  those  points  from  V.  In  the  case  of  a  PAR,  it  isremoved  by  dropping 
its  vertex  point  (k)  from  V.1K  where  condition  (1)5  is  checked  by  sensing  on  s2  s  sin2  A0k  <  to2. 
(See  (47).  page  28).  When  this  inequality  is  satisfied  a!  ak  we  say  points  (k  -  I).  (k).  (k  +  1)  are 
successive  colincar  (mints  (SCP).  Note  this  inequality  is  also  satisfied  if  (Af^l  <  to.  so  vertex  points 
of  such  angular  regions  are  also  dropped.*8 

D  Definitions  here  for  PAR  and  SDP  arc  slightly  changed  from  thoc  given  00  page  12,  to  account  for  the  finite 
precision  of  the  CDC-6700. 

*8Tho  value  of  P(H)  is  not  chaugcd  by  dropping  such  points. 
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After  a  point,  or  a  succession  of  points,  are  deleted  from  V,  it  is  “closed  up”  (CU).  This  means 
the  array  is  brought  together  so  that  no  gaps  occur  with  the  points  renumbered  in  order.  For 
example,  in  Figure  4,  4,  5.  6,  7,  8  would  be  dropped  from  V  by  SORT  III,  3,  and  2  or  4  would 
evaluate  P  over  the  CU  array  ( 1 , 2, 3, 9,  10,  11,  1 2),  which  we  again  call  V.  Thus,  3  must  not  only 
detect  when  (A)  or  (B)  holds,  but  it  must  also  delete  points  from  V  and  CU  the  array.  In  giving 
the  details  of  3,  we  make  use  of  its  flow  chart  on  page  43. 

At  each  stage  of  3,  an  angular  region  a k  is  under  examination.  It  is  made  up  of  a  vertex  point 

(k) ,  a  preceding  point  (k-  1),  and  a  following  point  (m)  (initially  m  =  k+  1).  Points  (k-  11,  (k), 
(m)  refer  to  their  order  as  listed  in  V,  where  V  may  no  longer  be  the  original  array  due  to  previous 
deletions. 

The  program  3  is  started  with  a,.  k  =  1,  m  =  2,  i.e.,  with  points  (N),  (1),  (2),  3-2.  If  (N)  and 

(l)  are  SOP,  then  (N)  is  dropped  from  V  by  setting  N  =  N  -  1,  3-8,  9.  This  is  repeated  until  (N) 
and  (1)  are  not  SDP.  Similarly,  (1)  and  (m)  are  tested.  If  they  are  SDP,  m  is  increased  by  one 
(m  =  m  +  1),  and  (1)  and  (m)  are  tested,  where  (m)  now  would  refer  to  (3)  of  V.  This  is  repeated 
until  ( 1 )  and  (m)  are  not  SDP,  3-15,  1 6.  The  array  V  is  then  reduced  by  deleting  the  proper  points 
and  then  CU  by  replacing  points  of  V  starting  at  (21  by  points  (m)  through  (N).  Then  N  is  replaced 
by  N  -  m  +  2,  3-17,  22.  The  value  of  N  now  refers  to  the  number  of  elements  in  the  updated  CU 
array  V.  Assuming,  at  this  point,  that  (A)  is  not  satisfied  fora,,  i.e.,  neither  (N)  ami  (1)  nor(l) 
and  (2)  are  SDP  in  V.  then  condition  (B)  is  checked,  3-18.  If  (N),  (1),  (m){=2)  are  SCP(s2  <to2), 
then  m  is  increased  by  one  until  (l)  and  (m)  are  not  SDP  and  (N),  (1),  <m)  are  not  SCP,  3-19, 
20,  1 8.  If  m  >  N.  3-19,  then  all  points  are  colinear.  N  is  set  to  2  and  3  exits  to  P-2.  If  2  <  m  <  N, 
3-23,  V  is  reduced  and  CU  by  replacing  elements  starting  at  ( 1)  with  elements  (m  -  1 )  through  IN), 
fhe  updated  V  will  now  contain  N  =  N  -  m  +  2  elements.  3-23,  27.  Control  is  returned  to  3-2. 

If  1  <  k  <  N,  and  if  tk)  and  tm)  are  SDP,  where  m  =»  k  +  1  initially,  then  a  new  angular  region 
is  considered  by  increasing  m  by  one  (m  -  m  +  1),  3-30  until  (k)  and  tm)  ate  not  SDP,  3-29.  Hie 
V  a  tray  is  then  reduced  and  CU  if  m  >  k  +  1 ,  by  replacing  elements  starting  at  (k  +  1)  by  elements 

(m)  through  (N).  The  updated  V  now  contains  N  3  N  -  m  +  1  +  k  elements,  3-28, 34.  Once  the  (A) 
condition  does  not  hold,  it  is  possible  to  check  if  M,  subtended  by  uk,  satisfies  II.  3*35,  If  (B) 
does  not  hold  and  uk  is  WD.  k  is  increased  by  one  (k  -  k  ♦  1)  and  the  procedure  for  1  <  k  <  N  is 
iterated.  3*39,  21. 


If  (B)  holds,  then  tk  -  1 ),  (k),  (m)  are  SCP  and  in  is  increased  by  one.  3*36.  The  value  of  m  is 
again  increased  by  one  if  (k)  and  fm)  are  SDP.  litis  is  continued  until  (k)  and  (in)  are  not  SDP  so 
that  (B)  can  be  tested  again.  3*36.  37,  35.  Eventually  (A)  and  (B)  arc  ruled  out,  if  in  >  k  +  1,  then 
V  is  reduced  and  CU  by  replacing  elements  starring  at  (k)  by  elements  (m  “  1)  through  iN).  Hie 
ululated  V  now  contains  N  *»  N  “  in  ♦  1  +  k  elements.  3*14.  (Note  if  in  *  k  ♦  1  at  3-39.  then  (A) 
and  (B)  do  not  hold  and  control  goes  to  3-21  to  look  at  (he  next  angular  region.)  If  at  3-36  m  >  N. 

then  k.  k  +  1.  ....  N  arc  colinear.  At  3-32  the  kih  point  is  replaced  by  the  N,h  point  and  N  is 

replaced  by  k.  The  updated  array  will  now  have  k  elements.  Control  is  passed  to  3-iO. 

Since  (k)  has  been  deleted,  (k  -  1)  and  the  new  (k)  clement  could  be  SDP.  If  they  are  ot. 
(hen  m  is  set  to  k  +  1  and  the  procedure  described  above  for  I  <  k  <  N  is  repealed  until  k  -  N. 

3*21.  30,  y.  If.  however  (k  -  I)  and  (k)  arc  SDP.  3-14.  then  k  is  reduced  by  one  (k  =  k  -  l). 
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3-38,  m  is  set  to  k  +  1, 3-33,  and  the  procedure  for  1  <  k  <  N  is  repeated  until  k  =  N,  3-21, 3-30,  or 
3-36.  In  the  event,  when  k  is  reduced,  that  it  takes  the  value  one,  the  entire  procedure  is  restarted 
with  che  updated  V  at  k  =  1,  m  =  2,  3-38,  3-2. 

1 

When  k  =  N,  the  Nth  angular  region,  specified  by  (N  -  1),  (N),(l)  with  respect  to  the  updated  V 
is  examined.  It  is  treated  in  much  the  same  way  as  a2.  The  details  may  be  gleaned  from  3-30,  36, 
31,32,21,10,11,5,12. 

A  final  possibility  exists  that  the  Nth  point  (N)  used  to  make  up  a  previously  WD  a j  is, 
subsequently,  deleted.  Therefore,  after  aN  has  been  accepted  as  WD,  ax  is  checked  again  to  assure 
that  (N),  (1),  (2)  are  not  SCP,  34,  6.  If  they  are  not,  then  an  exit  is  made  to  P-2  with  the  updated  V 
available  to  VALR-2  or  VALR-7  depending  on  the  value  of  ICV.  If,  however,  they  are  SCP,  then 
the  entire  procedure  starting  with  k  =  1 ,  m  =  3  must  be  carried  out  again  with  the  updated  V, 
3-7,  20.  This  takes  place  in  the  decomposition  of  the  element  in  Figure  30. 

Although  the  ideas,  and  general  description,  given  here  appear  straightforward,  their  imple¬ 
mentation  into  a  computer  program  that  handles  general  situations,  i.c.,  for  any  element  of 
{S}  using  3  with  4,  or  for  any  element  of  {11}  using  3  with  2,  requires  an  intricate  code  which  is 
reflected  in  the  flow  chart  of  SORT  III. 

If  3  is  used  with  3,  there  is  some  duplication  of  effort,  since  both, routines  check  for  SDP.  Of 
course,  2  can  be  used  alone,  as  mentioned  before,  for  any  II  in  {II},  but  it  may  be  more  efficient 
to  use  3  with  2  if  IT  contains  many  PAR(s),  since  an  erfc  function  computation  is  required  for  each 
of  them  when  2  is  used  alone,  which  does  not  occur  if  SORT  III  is  used  first. 

We  elaborate  the  discussion  of  3  by  processing  the  polygonal  elements  of  Figures  30  and  3 1  on 
pages  34-37.  The  element  in  Figure  30  is  in  {S}.  The  element  in  Figure  31  is  also  in  {S},  but 
computationally,  on  the  basis  of  the  discussion  on  pages  22, 23  (Figures  22  and  23),  both  30  and  3 1 
must  be  considered  as  self-intersecting.  We  shall  refer  to  both  of  them  as  S.  Their  processing 
involves  every  box  of  Flow  Chart  3. 

The  description  is  given  in  tabulated  form  on  pages  34-37.  The  first  column  contains  N.  the 
number  of  elements  in  V  at  certain  stages  of  the  processing.  The  second  column  lists  the  value  of  k, 
where  (k)  denotes  the  klh  node  or  point  of  V.  It  refers  to  the  nodal  point  (k)  which  with  (k  -  1 ) 
and  (m)  define  ak.  The  value  of  m  is  shown,  at  intermediate  stages,  in  column  3.  The  fourth  column 
displays  the  numbers  of  flow  chart  boxes  in  the  order  they  are  processed.  Column  five  shows  which 
points  are  deleted  from  V.  The  letters  preceding  the  dropped  points  are  helpful  to  establish  the 
updated  version  of  V  after  deletions,  For  example,  (a)  at  the  head  of  the  sixth  column  refers  to  the 
original  V  with  N  =  24  for  S  and  the  seventh  column,  headed  (b),  represents  the  CU  array  after 
points  (l),  (2),  (3),  shown  in  the  fifth  column  have  been  deleted  from  V  as  shown  under  (a).  Note, 
element  (4)  of  the  original  V  is  the  first  eler  ent  of  the  updated  CU  array  V.  in  column  (b).  which 
now  contains  21  elements,  i.e„  N  =  21  at  that  stage.  The  updated  V  at  each  stage,  where  one  or 
more  points  is  dropped,  is  shown  in  CU  form  by  columns  (b),  (c), ....  (q).  The  numbering  of  the 
elements  in  these  columns  retains  the  original  numbering  of  the  elements  in  V. 
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PROGRAM  3  FOR  S  FROM  FLOW  CHART 


BASED  ON  FIGURE  30 


N 
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m 
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PROGRAM  3  FOR  S  FROM  FLOW  CHART  (Continued) 
BASED  ON  FIGURE  30 
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PROGRAM  3  FOR  S  FROM  FLOW  CHART  (Continued) 
BASED  ON  FIGURE  30 
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PROGRAM  3  FOR  S  FROM  FLOW  CHART 
BASED  ON  FIGURE  31 
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As  mentioned  earlier,  there  is  no  need  to  discuss  4  (VALR-7)  extensively,  because  much  of  the 
coding  in  2  (VALR-2)  is  shared  by  4  (VALR-7).  Routine  4  yields  P(I1)  when  n  =  S  or  S  with  no 
SAR(s),  such  as  in  Figures  35  and  36.  Of  course,  pre-processing  S  by  3  to  remove  SAR(s)  allows  4 
to  be  used  with  the  reduced  element  which  will  then  be  in  {S}.  Routine  2  is  more  robust  than  4, 
because  it  can  find  P(II)  for  any  element  in  {II}.  Routine  4,  when  it  can  be  used,  is  preferred  to  2, 
because  it  may  be  more  efficient.  This  requires,  however,  that  the  user  must  know,  a  priori,  that  he 
has  an  element  in  {S}.  VALR-7  cannot  be  used  alone  for  an  element  in  {5}  with  SAR(s),  and 
cannot  be  used  at  all  for  an  SI  element.  If  it  is  used  alone  under  any  of  these  circumstances,  the 
value  for  P  will  most  likely  be  wrong. 

Some  specific  differences  between  4  and  2  are: 

(1)  4  uses  SMP-7  to  compute  A,  4-9.  The  sign  of  A  determines  whether  (24)  or  (26)  is  used 
to  find  P,  4-58,  59,  62.  The  quantity  A  is  computed  internally  in  2  as  a  by-product;  it 
has  no  specific  use  there,  2-22,  68. 

(2)  The  setting  of  the  error  indicator  IND  is  normally  set  to  zero.  If  IND  =  2  in  2  or  4,  then 
the  input  polygon  contains  a  PAR.  In  this  case  the  output  for  2  is  correct,  but  for  4, 
P  and  A  are  probably  wrong.  If  IND  =  1  in  2,  then  the  input  N  was  specified  as  one  and 
SDP  occurred.  If  IND  =  3  in  2  or  4,  then  the  input  N  <  1  or  N  =  2.  For  IND  =  1  or 
IND  =  3  the  output  is  meaningless. 

(3)  A  winding  number  W  is  not  computed  in  4  ((38)  and  (39)),  since  4  never  treats  an  SI 
element  alone.  This  has  the  advantage  that  if  A6  is  not  used  for  a  particular  ak,  say  if  R 
is  large  or  sin  AO  is  small,  a  call  to  the  arctangent  subroutine  can  be  bypassed,  4-20,  18, 
1 1, 12,  13,4,  5,  6,  7,  35,  38.  For  2,  all  A 0  must  be  computed  to  evaluate  W  which  is 
needed,  since  2  is  based  on  (32),  Thus  4  should  be  used  instead  of  2  for  simple  polygons, 
and  also  for  S  elements  without  SAR(s). 

Finally,  it  is  recalled  that  a  polygon  may  often  be  specified  by  either  of  two  numbering 
schemes  called  a  and  0-options  (see  page  24).  Generally  0  is  the  desired  option  for  computation, 
since  it  may  require  fewer  points  to  specify  the  given  polygon  (see  page  24),  However  if  the  user 
wishes  to  determine  beforehand  the  class  to  which  a  polygon  belongs,  it  should  always  be  numbered 
using  the  a-option.  See  pages  18-25  for  more  discussion. 


Definitions  and  Page  References  for  Flow  Chart  Quantities 

IBND  -  Integral  of  the  bivariate  normal  density  function,  page  2 
P  -  Value  of  P-function  for  a  polygon  or  an  angular  region,  page  1 

A  -  Value  of  A-function  for  a  polygon,  page  9,  Appendix  D 

ICV  -  Program  input  parameter,  page  25,  Appendix  F 
IND  -  Program  output  error  parameter,  page  26,  Appendix  F 
IOP  -  Program  input  accuracy  parameter,  page  25,  Appendix  F 
W  -  Winding  number,  page  1 8 

al>  a2>  a3>  R2/ 2  -  Program  parameters  which  depend  on  IOP  -  Appendix  E 

B  =  R2/2,  page  28 
gj  =  R/\/2  cos0j,  i  =  1,  2,  page  6 
ht  =  R/\/2  sin  0i(  i  =  1,2,  page  6 

G  =  2v^  (tl2  _hl)  ”  i  (g2h2  -giV,  Page  7 

u,  v,  w,  z  -  Defined  in  Flow  Charts  2  and  4 

am_l  =  Chebyshev  coefficients  for  erfc  (u),  page  6 

gj  =  7  X  10"14 ;  o  s  5  X  10“14  (used  in  SORT  I),  Appendix  E 

E(h)  s  erfc  (h),  page  5 

E(h)  s  erf (h),  page  28 

f l  =  Multiple  of  ±2ir,  page  20 

AO  =  tair1  (^/$)  page  5  (See  Footnote  P,  page  20) 

4>  3  VW  -  uz 
<f>  3  UW  +  V2 

s  =  sin  AO  =  2i/,/D1D2,  page  28 
D,  =  (2(wJ  +  ?2)),/2 
D,  »  ( 2(u?  +vQ),/2 
c  =  (R2/2)cosA0  =  g,g2  h|h2 

Program  Identification  Number: 

P-2  1  P-7  5 

VALR-2  2  SORT  I  6 

SORT  III  3  SORT  II  7 

VALR-7  4  SMP-7  8 
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IW  CHART  2  (Continued) 
VALR-2  (Continued) 


FLOW  CHART  3 


Output:  P,  A,  IND 


exp  f — f  B  +  In  jt|  ) 


FLOW  CHART  4  (Continued) 


VI.  NUMERICAL  RESULTS 
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In  this  section,  a  variety  of  polygon  elements  are  shown  through  Figures  30-58.  Each  element 
has  its  vertices  and  relevant  edge  intersections  numbered  sequentially,  from  1  to  N,  in  the  order  it  is 
generated.19  The  numerals  are  not  always  placed  optimally  for  viewing;  nevertheless  the  xy- 
coordinates  to  which  a  numeral  corresponds  can  always  be  decided.  A  numeral  (node)  is  generally, 
but  not  always,  located  below  and  slightly  to  the  left  01  the  xy-point  to  which  it  belongs.  The 
xy-coordinate  values  are  always  rational  numbers  that  can  be  read  from  the  figures  using  a  tuler 
graduated  in  tenths.  Each  figure  lists  the  following  information: 

ICV,  P,  A,  Classification  ( _ ),  N. 

All  the  computations  were  performed  using  P-2  as  the  master  routine.  If  there  is  a  most 
efficient  way  to  compute  P,  ICV  is  assigned  one  value.  If  there  are  two  ways  which  may  be  equally 
good,  or  if  it  is  difficult  to  decide  which  is  better,  then  ICV  is  assigned  two  values.  For  example. 
Figure  37  shows  a  simple  polygon,  so  ICV  =  0.  A  glance  at  the  P-2  flow  chart,  page  40,  confirms 
that  P-2  calls  VALR-7  to  compute  P(S).  In  Figure  32,  an  SI  element  is  shown  with  a  PAR  as  well. 
We  have  ICV  =  1,-1,  which,  from  P-2,  calls  first  VALR-2  to  find  P  and  subsequently  VALR-2 
preceded  by  SORT  III  to  obtain  the  same  result.  The  rounded  value  of  P  given  in  each  figure  is 
correct  to  the  number  of  digits  shown.  The  value  of  A  is  given  next.  The  classification  of  an 
element  (according  to  the  T-construction,  page  15)  follows  and  is  designated  by  S,  S,  or  SI.  The 
two  blanks,  in  parentheses,  following  the  classification,  as  noted  above,  are  used  to  denote  the 
computed  winding  number  W.  It  is  only  listed  if  VALR-2  is  called.  Thus,  for  Figure  32  two  winding 
numbers  are  listed  since  ICV  =1,-1  which  both  use  VALR-2.  In  the  first  case  the  arctangent  sub¬ 
routine  yields  ir  for  the  angular  measure  of  the  PAR  at  (11)  instead  of -it  according  to  the 
T-construction.  This  results  in  a  computed  winding  number  of  one  instead  of  mo.  Sec  pages  22 
and  23,  Note  there  is  also  a  PAR  at  ( 1 9);  however  in  this  case  W  is  not  affected  since  it  is  its  measure 
according  to  the  T-construction.  In  Figure  3?,  W  =  1.  since  the  element  is  PC),  but  it  is  not  listed 
because  VAI.R-2  was  not  used  to  compute  P  for  this  element.  Finally.  N  is  listed  which  refers  to 
the  number  of  points  used  to  define  the  configuration  as  shown, 

By  our  T  construction,  page  15.  an  element  of  (S)  has  a  winding  number  W  of  •.  I  However, 
because  of  the  range  of  the  arctangent  routine,  page  22.  this  need  not  be  the  ease  computationally 
as  for  example  in  Figures  30  and  31.  see  page  33  also,  lire  element  of  Figure  30(31)  is  in  (S'<  ac¬ 
cording  to  the  T-construction,  but  must  be  considered  J>1  for  computations.  Thus  P  is  computed 
using  VALR-2  with  a  computed  W  of  6(2).  using  ICV  *  l.  Then  P  is  computed  again  using  VAlR-7 
preceded  by  SORT  III.  (ICV  =  -2). 

In  Figure  34,  we  have  the  case  ot  a  simple  polygon  in  the  form  of  a  triple  spiral.  Another 
simple  polygon  is  shown  in  Figure  37.  Figures  30,  31. 35.  3b.  38.  39,  42.  43,44.  45.48.  49.  51. 
55.  5 6  contain  elements  in  (S’.  Die  remaining  figures:  32.  33. 40.  41. 4b.  47.  50.  52.  53.  54.  57. 
58  display  SI  elements 


*9lt  should  be  understood  that  an  additional  node  (S'  ♦  I)  is  located  ai  MS  ‘  1 ). 
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It  is  our  objective  in  presenting  these  figures,  that  there  is  enough  variety  to  resolve  for  the 
reader  any  remaining  uncertainties  regarding  the  robustness  of  (P-B),  how  a  polygon  is  specified 
and  classified,  and  how  a  winding  number  is  determined.  Finally,  for  completeness  and  as  a  further 
clarification  of  the  role  of  exterior  angular  regions,  we  tabulate  on  page  48  a  detailed  listing  of  the 
values  of  P(ak),  k  =  1,  2,  N,  that  are  needed  to  compute  P(fl)  for  the  element  displayed  in 
Figure  33. 

This  polygon  is  interesting  in  its  own  right,  since  P(I1)  represents,  here,  the  probability  for 
an  event,  governed  by  a  bivariate  normal  distribution,  occuning  in  Sj  and/or  S2,  where  S]  = 
(1,  2,  3,  4,  5,  6,  7),  S2  =  (8,  9,  10,  1 1,  12).  From  probability  theory,  we  can  write 

(54)  p(ri)  =  pcs!  us2)  =  P(Sj)  +  P(S2)  -  |P(S,  ns2)i, 

where  U,  Pi  denote  union  and  intersection  symbols  for  sets  and  Sj  Pi  S2  =(8,  13,  14,  15,  16,  17,  18) 
which  is  NO.  The  values  of  P(Sj ),  P(S2),  P(Sj  0  S2)  are  given  in  Appendix  A,  where  P(H)  is  found 
by  decomposing  n  with  SORT  i,  which  is  based  on  (P-A).  See  page  (A- 14). 

The  tabulation,  page  48,  lists  in  the  first  column  the  value  of  k  for  the  kth  node  of  IT  The 
second  and  third  columns  list  the  x  and  y  coordinates,  respectively,  for  each  node  in  the  order  n  is 
generated  (the  numbering  used  is  under  the  /3-option,  see  page  24).  The  fourth  column  lists  the 
value  of  P(ak),  for  each  f.  obtained  from  VALR-2  with  1CV  =  1  and  IOP  =  3  in  P-2.  The  corre¬ 
sponding  angular  measures  A0k  for  each  ak  are  given  in  t’.ie  fifth  column.  In  the  sixth  and  seventh 
columns  P(ak)  and  A0k  are  listed  for  ICV  =  -1,  i.e.,  n  is  treated  by  SORT  III  first,  and  then  values 
in  the  sixth  and  seventh  columns  are  computed,  with  IOP  =  3,  from  VALR-2.  (Note,  IOP  =  3 
implies  an  accuracy  of  9-decimal  digits  for  each  P(ak ).)  The  reduced  polygon  as  a  result  of  SORT  III 
treating  II  is  given  by  (1,  2,  3,  4.  5,  6,  7,  8.  9,  10,  14,  15, 16,  17, 18, 19).  SORT  III  has  deleted  (12) 
since  (11),  (12),  (13)  are  SCP.  Then  it  drops  (13),  because  (11)  and  (13)  are  now  SDP.  Then  it 
removes  (11)  because  (10),  (11),  (14)  are  now  SCP.  The  columns  4,  5,  6,  7  are  summed  at  the 
bottom  of  the  tabulation.  Note  that  W  =  2  for  both  situations.  This  is  not  always  the  case.  The 
primary  circuits  (see  page  18)  can  be  gleaned  from  the  figure.  For  ICV  =  l,  the  circuits  are  given  by 
CpOI,)  -  (16,4,5,6,7,8,  9,  16)  with  W,  =  1,  C,,(II2)  =  (1,2,3,  16,  iO.  11,12,  13,14,  15,  16,  17, 
18,  19)20  with  W2  =  1.  Hence  W  =  Wt  +  W2  =  2.  For  ICV  =  - 1,  the  primary  circuits  are  given  by 
CP(n,)  =  (16,4,5,6,7,8,9,  16)  with  W,  -  1.  CP(I12)  =  (1.  2,3,  16, 10, 14.  15,  16,  17,  18, 19).  with 
W2  -  1.  He.'ce  again  W  =  2.  A  7r-angular  region,  PAR,  occurs  at  k  =  12  initially.  Hence  with  ICV  =  1 
an  erfc  calculation  is  required  at  k  =  12,  namely  1/2  erfc  (h2)  *  1/2  erfc  (~1)  =  .92135  03965,  If 
ICV  =  -1,  then  point  (12)  is  deleted  and  the  erfc  computation,  at  the  expense  of  using  SORT  HI. 
is  not  necessary. 

The  time  of  computation  per  angular  region  is  given  in  Appendix  E. 

All  of  the  numerical  results  in  this  report,  as  well  as  many  that  are  not  given,  were  checked 
by  an  independent  procedure.  It  consists  of  decomposing  II,  regardless  of  its  class,  into  a 
set  of  triangles  {Aj }.  The  triangle  Aj  has  vertices  (1 ).  (j),  (j  +  l ).  With  j  =  2.  3 . N  -  I  we  have 


20The  order  of  the  nodes  appears  unusual  because  of  the  use  of  the  0-option  numbering  scheme. 


P(IT)  =  Sjlj1  P(Aj).  The  proof  of  this  result  follows  the  lines  of  proof  given  for  A  in  Appendix  D; 
it  is  not  given  here. 

The  value  of  P(Aj)  is  computed  from  a  routine  we  developed  which  uses  Drezner’s  scheme 
[2,  page  18]  for  evaluating  P  ewer  an  angular  region.  His  method  is  much  slower  than  ours  but 
gives  very  good  accuracy.  It  is  described  in  [2] . 

In  this  checkout  program,  which  is  listed  in  Appendix  G;  N-2  triangles  A:  are  obtained  for 
each  II,  and  consequently  P  is  required  for  3(N-2)  angular  regions. 


TABULATION  OF  RESULTS  FOR  FIGURE  33  USING  P-2,  (e  =  5  X  1(T10) 


"  "'"'l 

k 

X 

— 

y 

ICV=1,  P(ak) 

ICV=1,  Aflk 

ICV  =  -1,  P(ak) 

ICV  =  -1,  A0k 

1 

-3 

0 

1.6803  81909  (-2) 

3tt/4 

1.6803  81909  (-2) 

37r/4 

2 

0 

-3 

7.8697  69659  (-3) 

1.4288  99272(0) 

7.869/  69659  (-3) 

1.4288  99272(0) 

3 

4 

0 

4.9999  79129  (-1) 

2.4980  91545  (0) 

4.9999  79129  (-1) 

2.4980  91545(0) 

4 

0 

0 

-3/8 

-3tt/4 

-3/8 

-3rr/4 

5 

2 

2 

3.0600  67674  (-3) 

1 .8925  46881  (0) 

3.0600  67674  (-3) 

1.8925  46881  (0) 

6 

0 

3 

1.6551  276 10  (-2) 

1.2490  45772(0) 

1.6551  27610  (-2) 

1.2490  45772(0) 

7 

-3 

0 

4.9985  63923  (-1) 

3rr/4 

4.9985  63923  (-1) 

3tt/4 

8 

-2 

0 

-3.1610  42924  (-1) 

--4.6364  76090  (-1) 

-3.1610  42924  (-1) 

-4.6364  76090  (-1) 

9 

0 

"1 

9.5914  28393  (-2) 

9.2729  52180  (-1) 

9.5914  28393  (-2) 

9.2729  52180  (-1) 

10 

4 

1 

1 

1.5865  448 1 3  (-1) 

2.6779  45045  (0) 

1.5865  44813  (-1) 

2.6779  45045  (0) 

It 

-1 

2.6739  05696  (-2) 

tt/4 

12 

0 

9.2135  03965  (-1) 

IT 

13 

-I 

1 

-1.0674  47074  (-1) 

-it/ 4 

14 

i 

1 

-4.8741  42552  (-1) 

-3*/4 

3.5393  04909  (-1) 

n/4 

15 

0 

0 

3/8 

3t*/4 

3/8 

3jt/4 

I  ft 

2 

0 

-1.8389  57076  (-1) 

-2.6779  45045  (0) 

-1.8389  57076  (-1) 

-2.6779  45045(0) 

17 

0 

-1 

-9.5914  28393  (’■2) 

-9.2729  5.180  (-1) 

-9.5914  28393  (-2) 

-9.2729  52 180  (-1) 

18 

"2 

0 

1.6509  77199  (-3) 

4.6364  76090  (-1) 

1.6509  77199  (-3) 

4.6364  76090  (-1) 

P(  11)  -  2  -•  21  P{<ik )  •  0.94 1 62  48 1 30 

(£A0k/2rr)  =  W»2 

P(»)  “0.94162  48129 

(21  Aflk/2ir]  «  W  =  2 
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APPENDIX  A 

NORMAL  PROBABILITY  OVER  ARBITRARY  POLYGONS  BY  (P- A) 


In  this  appendix,  we  describe  Procedure  (P-A)  for  computing  P(II),  n  in  {II}.  It  differs 
significantly  from  (P~B)  which  is  discussed  in  Sections  IV,  V  and  VI.  In  (P  —  B),  the  concept  of  a 
winding  number  was  introduced.  A  program  package,  called  Program  B,  was  developed  made  up  of 
five  subprograms:  P-2,  VALR-2,  SORT  III,  VALR-7,  and  SMP-7.  The  (P- A)  procedure,  like  (P-B), 
was  developed  to  treat  SI  elements  or  S  elements  with  SAR(s).  It  decomposes  such  an  element  into 
a  set  of  disjoint  elements  in  {S }  or  { S }  depending  upon  whether  the  (3  or  a-option  is  used  to  specify 
n  (See  page  24).  If  the  decomposed  elements  are  in  {S}  they  are  treated  by  removing  any  PAR(s). 
The  resulting  decomposed  elements  are  then  in  {S},  and  P  for  such  elements  is  computed  by  using 
VALR-7.  A  program  package  is  also  described  in  this  appendix,  using  (P- A),  which  is  comprised  of 
six  subprograms:  P-7,  VALR-7,  SORT  I,  SORT  II,  SORT  III  and  SMP-7.  This  program  package  is 
called  Program  A. 

Procedure  (P-  A)  has  merit,  because  its  decomposition  of  H  into  disjoint  elements  in  { S }  (or 
{S})  allows  the  analyst  to  gather  a  more  detailed  picture  of  the  type  of  region  n  represents.  More¬ 
over  since  it  permits  the  decomposition  to  be  carried  out  in  a  backward  (from  (N)  to  (1))  as  well  as 
forward  direction,  it  gives  an  additional  means  of  checking  final  results  (The  decomposition  is  not 
necessarily  the  same  in  the  forward  and  backward  directions).  Nevertheless  (P-B)  is  deemed  the 
better  overall  procedure.  Some  summarizing  remarks  comparing  (P-A)  and  (P-B)  arc  given  on 
page  A- 15. 

For  computational  efficiency  with  Program  A,  the  pre-processing  of  II  must  be  kept  to  a 
minimum  by  specifying  beforehand  the  smallest  class  to  which  11  belongs.  However  if  fl  is  errone¬ 
ously  assigned  to  a  class  to  which  it  does  not  belong,  then,  because  VALR-7  cannot  treat  SI  elements 
nor  elements  in  {S}  with  SAR(s)  a  wrong  result  for  P(II)  is  likely.  Recall,  on  the  other  hand,  that 
VALR-2  of  Program  B  can  treat  any  polygon  with  the  same  computational  efficiency,  although 
small  improvements  in  efficiency  can  sometimes  be  achieved  by  pre-processing  II  with  SORT  III 
and  by  using  VALR-7  for  VALR-2  as  indicated  in  P-2  (see  page  40). 

We  assume  throughout  this  appendix  that  (P~A)  or  Program  A  are  under  discussion  unless 
stated  otherwise. 

If  II  is  in  {S},  then  by  computing  P(ak)  for  each  «k  of  S,  P(S)  is  obtained  by  using  (24)  if 
A  >  0,  and  (26)  if  A  <  0.  In  the  first  case  n  is  PO  and  in  the  second  II  is  NO. 

If  II  is  in  { S } ,  then  SAR(s)  can  occur  (see  pages  12,  31).  The  SDP  and  SCT  are  treated  by 
removing  appropriate  points  from  the  V-array  which  specifies  II.  The  reduced  polygon  is  in  { S) , 
and  P(S)  (=P(S))  is  obtained  as  explained  in  the  previous  paragraph. 

If  n  is  SI,  then  II  is  decomposed  into  a  set  of  disjoint  elements,  often  referred  to  as  isolated 
elements ,  in  {S}  or  {S}  depending  on  whether  11  has  been  specified  by  the  a  or  (3  numbering  scheme. 
If  II  is  decomposed  into  (S1  •••  Sn)  then  P(II)  =  2"  PfS1)  and  if  the  decomposition  results  in 
(Sl  •••  Sn),  then  each  S1  must  be  processed  first,  as  explained  in  the  preceding  paragraph,  before 
P(S!)  can  be  computed  such  that  P(II)  =  £('  PIS1). 
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The  a  and  /3-options  come  into  play  if  II  is  in  {S>  or  {11}  (see  page  14).  Generally  the  /3-option 
is  the  more  efficient,  since  it  will  often  require  the  treatment  of  fewer  angular  regions  (one  must 
be  careful,  however,  that  it  specifies  II  correctly).  In  Figure  A-l,1  16  angular  regions  occur  with 
the  a-option.  For  the  same  II  using  the  /3-option,  as  shown  in  Figure  A-2.  only  ten  regions  occur, 
and  two  of  those  are  of  no  consequence  since  A8  =  0  for  them. 


Figure  A-l.  Polygon 
with  a-Option 


Figure  A-2.  Polygon 
with  /3-Option 


The  two  figures  above  also  show  that  to  determine  whether  II  is  SI  every  vertex,  meeting  or 
intersection  of  two  edges  must  be  numbered  each  time  it  occurs  in  the  order  it  occurs  (see  page  14). 
Subsequently,  the  polygon  can  be  numbered  with  the  /3-option,  if  the  opportunity  exists  to  do  so, 
for  the  actual  computation  of  P(H).  Note  that  although  n  is  SI  in  Figures  A-l  and  A-2,  this  would 
not  be  concluded  from  our  T-characterization  (page  15)  by  examining  the  numbered  points  in 
Figure  A-2. 

If  II  is  SI  and  numbered  under  the  a  or  0-options,  then  by  (P-A)  II  is  decomposed  in  the 
following  way: 

Starting  at  node  (1),  we  look  for  the  first  MN  (see  page  14)  that  is  met  for  the  second  time, 
say  MN(k)  is  met  for  the  second  time,  say  at  node  (k  +  m),  1  ^k«$N"l,2<k  +  m<N  +  1, 
(k  +  m  =  N  +  1  means  (k  +  m)  =  ( 1 )).  In  Figure  A-3,  this  occurs  at  MN(3),  since  this  point  is  first 
encountered  by  node  (3)  and  for  the  second  time  by  node  (7).  The  same  property  also  holds  for 
Figure  A-4  at  MN(3),  which  is  met  for  the  second  time  at  node  (6). 

When  this  situation  occurs,  there  exists  two  possibilities,2  [4,  p.  16] : 

(a)  Edges  k  +  T  and  F+  m,  m  >  2,  with  k  +  1  originating  at  MN(k)  and  k  +  ni  terminating 
there,  have  mom  than  one  point  in  common,  (Under  a-option  m  =  2). 


1  Figures  A-l  anti  A-2  are  the  same  as  Figures  25  and  26. 

2  Actually,  a  third  possibility  oxists,  namely  SDP.  For  each  sot  oi'SDP,  as  soon  as  it  is  detected,  one  of  the  points  is 
removed. 
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Figure  A-3.  An  SI  Polygon  with  IBND 
Required  Over  Shaded  Areas,  N  =  1 9 


Figure  A-4.  An  SI  Polygon  with  Line 
Segments.  IBN  D  Required  Over  Shaded 
Areas,  N  =  22. 


(b)  Nodes  (k),  (k  +  l) . (k  +  m),  m  >  2,  specify  a  polygon  in  {S}  (Under  adoption  the 

polygon  is  in  { S> . )  with  a  clearly  defined  orientation  (see  page  9).  Here  also  k'V’l 
originates  at  MN(k)  and  k  +  m  ends  there. 

in  case  (a),  a  line  segment  of  II  exists  where  edges  k  +  1  and  k  +  m  completely  overlap.  The  con¬ 
figuration  of  Figure  A-4  contains  examples  of  such  line  segments.  They  will  be  identified  below 
where  II  in  that  figure  is  decomposed. 

In  case  (b),  an  element  of  (S>  or  is  obtained  with  M  =  m  nodes.  The  function  P  is  com¬ 
puted  for  this  element  (using  VALR-7).  The  nodes  (k)  to  (k  +  m  -  1)  are  then  deleted  from  the 
original  V-array  specifying  II,  and  the  decomposition  continues  starting  at  (k  +  ml,  which  is  now 
the  k,!l  element  of  the  updated  V-array.  Since  II  has  only  N  nodes,  this  will  end  after  a  finite 
number  of  such  steps.  P(ll)  is  computed  by  adding  up  the  positive  contributions  from  PO 
isolated  polygons  and  the  negative  contributions  from  the  isolated  NO  polygons. 

A  proof  that  the  above  decomposition  can  always  be  carried  out  is  essentially  given  in  Knopp, 
[4,  page  151 .  His  proof,  which  requires  minor  changes  for  our  use,  is  constructive.  We  have  used  it 
as  a  guide  in  the  decomposition  procedure  just  described. 
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In  detailing  the  decomposition  of  the  polygons  in  Figures  A-3  and  A-4,  the  isolated  simple 
polygons  are  superscripted  in  the  order  they  are  isolated,  i.e.,  S1,  S2.  ....  Sn.  They  are  identified, 
as  usual,  by  their  nodes.  We  also  give  their  orientation.  They  are  both  specified  by  the  a-option. 

Figure  A-3 

S!:  (3,4,  5,  6,  7)  PO  S3:  (12,  1 3,  14,  15)  NO 

S2:  (2,  7,  8,  9,  10,  11)  NO  S4:  ( 1, 1 1,  15,  16,  17,  18, 19, 20)  PO 

Note:  Case  (a)  does  not  occur  here. 


S1:  (3,4,  5,6)  PO 
S2:  (6,  7,8,9)  NO 
S3:  (11,  12,  13,  14)  PO 
S4:  (14,  15,  16,  17)  NO 


Figure  A-4 

Line  Segment:  (10,  17,  18)  Case  (a)3 
Ss:  (18,  19,20,21)  PO 
Line  Segment:  (9,21,22)  Case(a) 
S6:  (1,2,22,23)  PO 


An  automatic  formal  procedure  for  decomposing  an  SI  polygon  n  is  carried  out  by  listing  the 
integers  corresponding  to  its  ordered  set  of  nodes.  In  general,  after  S'  is  found,  i  n,  all  the  integers 
corresponding  to  the  nodes  of  S*,  except  the  last,  are  dropped  from  the  initial  list  V.  However  if 
S1  contains  node  ( 1 )  then  one  is  retained,  rather  than  the  integer  corresponding  to  the  last  node  of 
S'.  For  example,  for  Figure  A-4,  we  would  have  after  deleting  S1 : 

V  =  1,2,  6,  7,  8,  9,  10.  11.  12,  13,  14,  15.  16,  17.  18,  19,20,21,  22.23. 

Starting  at  6,  S2  is  found,  and  the  above  list  is  reduced  to 

V  =  1,  2.9,  10,  11.  12,  13,  14.  15.  16.  17,  18,19.20,21,22,23. 

After  S3  and  S4  are  found,  we  have 

(•)  1,  2,9,  10,  17,  18,  19,20,21.22,23. 

At  this  stage,  we  get  lino  segment  (10,  17,  18)  and  the  set  (•)  above  contracts  to 

V  =  1,  2.9,  18.  19.  20,21.  22,  23. 

After  S5  is  found,  we  have 


V 


1.  2,  9.  21.  22.  23. 


The  removal  >f  the  line  segment  (9,  21,  22),  leaves  us  with  (1.2, 
the  decomposition. 


22,  23)  which  is  S6  and  concludes 


Figures  A-5  and  A-6  contain  the  same  polygon  with  the  a  and  0-options,  respectively.  From 
the  details  of  the  decompositions  given  below,  it  will  be  clear  PAR(s)  under  the  a-option  are 
removed  during  the  decomposition,  but  they  can  be  retained  under  the  0-option  as  this  example 


3Note  that  the  decomposition  isolates  SAR(s),  if  the  ot-oplion  is  used. 
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Figure  A-5.  SI  Polygon,  a-Option,  N  =  13  Figure  A-6.  SI  Polygon, /3-Option,  N  =  1 1 

shows.  Thus,  with  the  /3-option  an  additional  program  SORT  II,  is  needed  to  eliminate  PAR(s). 
However  it  actually  does  a  little  more  by  eliminating  SCP(s)  (see  p.  31).  This  program  works  in 
the  same  way  as  SORT  III  only  it  need  not  cheek  for  SDP,4  (see  page  31). 

For  Figure  A-5,  the  decomposition  by  SORT  l  gives: 

S!  =  (2,  3, 4,  5)  PO 

Line  Segment  (l,  5, 6)  (Fliminated  by  SORT  1  since  a  PAR  requires  p  =  0) 

S'  =  (9.  10,  11, 12)  PO 

Line  Segment  (8,  12,  13)  (Eliminated  by  SORT  1  since  a  PAR  requires  p  =  0) 

SJ  =(1.7,  13.  14)  NO. 

With  Figure  A-6.  the  decomposition  by  SORT  1  yields: 

51  =  (2, 3, 4. 5)  PO 

52  =  (8,9.  10. 11)  PO 

SJ  =  (1,5.  6.  7,  11.  12)  NO 

where  points  (5)  and  ( 1 1 )  are  eliminated  from  SJ  by  SOR  T  li. 

Note  that  in  Figure  A-h,  it  was  necessary  to  include  the  point  (II)  coinciding  with (8) or (7). 
otherwise  the  decomposition  would  have  left  the  SI  polygon  (I.  5,  6.  7,  8.  9,  10.  12)  with  no 
coinciding  points.  Hence,  after  SORT  II,  which  would  remove  (5),  the  resulting  element 
d,  6,  7.  8,  V.  10,  12)  is  SI  and  VAI.R-7  applied  to  it  would  yield  a  wrong  result.  This  example 
serves  to  emphasize,  with  the  (P-A)  procedure,  the  care  that  must  be  taken  in  using  the  /3-option 
to  specify  11.  The  (P-B)  procedure  would  have  no  difficulty  in  this  situation  since  VALR-2  can 
handle  SI  polygons  directly. 

We  proceed  with  a  description  of  the  computer  program  package  based  on  (P-A),  i.e.,  Pro¬ 
gram  A.  Recall  that  for  (P-B),  Program  B  is  composed  of  P-2,  VALR-2,  SORT  HI.  VALR-7, 


4It  should  be  evident  that  SDP  arc  always  delected  and  removed  by  SORT  l  by  testing  if  M  <2  (sec  box  6  of  the 
Flow  Chart  6). 
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SMP-7.  For  easy  reference,  they  were  referred  to  as  subprograms  1,  2,  3,  4  and  8  respectively. 
For  Program  A,  the  program  package  consists  of  P-7,  VALR-7,  SORT  I,  SORT  II,  SORT  III,  and 
SMP-7.  For  easy  reference,  we  number  them  accordingly;  5  P-7,  4  ■**  VALR-7,  6  SORT  I, 

7  **  SORT  II,  3  SORT  III,  8  SMP-7.  All  of  these  subprograms  are  in  subroutine  format.  Pro¬ 
gram  P-7  serves  as  a  master  routine,  VALR-7  is  much  like  VALR-2,  but  it  can  only  compute  P  for  a 
single  angular  region,  or  polygons  in  {S},  provided  they  contain  no  SAR(s).  SORT  I  decomposes  n 
into  a  set  of  disjoint  elements  in  {S}  or  {S}  depending  on  whether  it  is  numbered  with  the  a  or 
0-option.  It  is  primarily  used  if  n  is  SI.  SORT  II  is  used  on  the  disjoint  elements  in  {S},  obtained 
from  SORT  I,  to  remove  SCP(s).  SORT  III  was  taken  up  in  Section  V.  It  is  used  to  delete  SC'P  and 
SDP  from  II,  the  original  polygon,  when  II  is  in  { S> .  SMP-7  is  used  to  compute  the  function  A  as 
given  in  (22);  where  1 A 1  is  taken  as  the  area  of  II;  the  sign  of  A  is  used  in  VALR-7  to  determine 
the  orientation  of  II  when  II  is  in  (S)  or  {§}. 

The  flow  charts  for  3  and  4  are  given  at  end  of  Section  V,  pages  43-45,  since  they  also  make 
up  part  of  Program  B.  Flow  charts  for  5,  6,  7  are  given  at  the  end  of  this  appendix,  pages  (A- 17- 
A-19).  No  tlow  chart  is  given  for  SMP-7.  Fortran  IV  listings  of  all  the  programs  are  given  in 
Appendix  F. 

Program  5  (see  Flow  Chart  5)  uses  as  input  x,  y,  N,  ICV  and  10P.  These  notations  have  all 
been  used  previously  in  Section  V.  The  various  values  for  ICV  have  slightly  different  meaning  here. 
If  ICV  =  0,  P(S)  or  P(S)  is  wanted  where  S  has  no  SAR(s)  such  as  in  Figure  35.  If  ICV  =  1,  then 
P(S)  is  wanted,  where  3  is  used  before  4  to  remove  SCP  and  SDP.  If  N  =  1 ,  P  for  an  angular  region 
is  wanted.  If  ICV  =  42  or  ±3,  P  for  an  arbitrary  polygon  is  wanted.  IOP  specifies  the  accuracy 
desired;  it  can  be  assigned  the  values  1,  2,  or  3  to  yield  approximately  3,  6.  or  9 -decimal-digit 
accuracy,  respectively,  for  P  of  each  angular  region. 

If  IlCVi  *  2,  it  is  assumed  that  the  op-option  has  been  used  to  specify  an  element  in  (Ill. 
If  IK'V'i  =  131,  it  b  assumed  the  0-option  has  been  used.  In  the  first  ease  the  isolated  elements  are 
in  (S)  and  in  the  second  case  they  may  be  in  (S).  If  ICV  =  2  or  3,  the  processing  of  11.  by  SORT  I, 
begins  at  (1)  and  progresses  sequentially  through  nodes  (2),  (3),  ....  (N).  If  ICV  *  -2  or  -3,  then 
ll  is  processed  by  SORT  I  in  reverse  order  starting  at  (N)  and  progressing  sequentially  through 
(N-  1UN-2) . (1). 

The  parameter  IND  is  discussed  below. 

Wc  now  consider  6  in  more  detail  by  using  its  fiow  cliart,  page  A- IS.  Two  points  of  (vt  e 
(xr  yr)l.  r*  1.  2 . K,  are  said  to  coincide  or  arc  duplicates  if 

(A-l )  |Xj-xk|  <  o.  lyj-ykl  <  o,  l  <  k  <  i  <  K,  o  a  5(-!4>. 

Program  6  is  started  by  setting  K  s  N  and  then  by  sensing  if  v,  and  vN  of  the  V-anay,  which  spec¬ 
ifics  II.  coincide .  6-3.  If  they  do,  then  v,  replaces  vN  in  the  V-array.  If  they  do  not  coincide,  then 
v,  is  added  to  V  as  vN+1  and  K  =  N  +  1.  Before  proceeding  with  the  decomposition  of  11, 6  deter¬ 
mines  whether  V  is  to  be  processed  in  increasing  order  of  its  elements  or  in  decreasing  order.  6-2. 
The  resulting  decompositions  are  not  necessarily  identical,  i.e.,  they  may  not  isolate  the  same  set  of 
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polygons.  Figure  36  contains  an  example.  The  two  decompositions  for  that  example  are  given 
near  the  end  of  this  appendix,  page  A-16.  Of  course,  the  result  for  P(I1)  must  be  independent  of 
which  decomposition  is  used. 

The  procedure  used  by  SORT  I  is  an  N2-process,  whereas  SORT  III  is  an  N-process. 

We  focus  our  attention  on  the  forward  decomposition,  (ICV  >  0),  64,  rather  than  the  reverse 
procedure,  (ICV  <  0),  6-5. 

The  array  V  of  data  points  is  searched  for  a  point  vk,  1  <  k  <  i,  which  coincides  with  v,, 
starting  with  i  =  2,  i.e.,  where  vs  and  vk  satisfy  (A-l).  If  vs  and  vk  coincide,  1  <  k  <  i,  then  set 
1ST  =  k  and  IEN  =  i,  with  2  -*  i  -*■  L  if  k  ■  1,  otherwise  k  -*•  i  L,  64.  In  6-6,  M  is  set  to  NUMI; 
the  inequality  IEN  -  1ST  =  NUMI  <  2  is  tested.  If  NUMI  =  1  or  NUMI  =  2,  then  we  have  isolated 
either  a  set  of  SDP  or  a  set  of  SCP,  respectively.  In  either  case,  such  elements  do  not  contribute  to 
P.  Hence  we  set  p5  =  0,  and  a  call  to  VALR-7,  6-9,  can  be  averted.  If  the  inequality  is  not  satisfied, 
then  an  element  of  (S)  or  {§}  has  been  isolated.  In  order  to  determine  whether  it  is  in  {S},  i.e., 
it  tire  ooption  has  been  specified,  a  sensing  on  ICV  is  carried  out  at  6-10.  If  the  answer  is  yes  at  this 
box,  then  the  isolated  element  is  assumed  to  be  in  {S},  (a-option),  and  SORT  II,  6-13,  is  not  ealled. 
It  the  answer  isjio,  then  SORT  II  will  be  called,  6-13,  since  it  is  assumed  in  this  case  that  the  isolated 
element  is  in  (SK  (0-option).  In  6-13,  the  inequality  NUMI  <  2  is  checked  again,  after  SORT  II 
Iras  been  used.  It  could  happen  that  after  deletions  by  SORT  II,  the  isolated  element  S  retains  no 
more  than  3  points,  so  that  the  inequality  NUMI  <  2  would  be  satisfied.  Then  p  =  0.  and  VALR-7, 
6-9,  is  bypassed;  the  program  proceeds  directly  to  64.  If  the  inequalities  of  6-6  and/or  6-13  are 
not  satisfied,  then  VAl  R-7  is  called  to  compute  P  and  A  for  the  isolated  clement,  which  we  denote 
here  by  p  and  a.  respectively. 

Following  tire  computation  of  p  and  a,  a  query  is  made  at  6-8.  Is  IEN  =  K?  If  the  answer  is 

no.  H  requires  further  processing,  which  is  carried  out  after  replacing  elements  L . K  -  M  of  V 

by  elements  (L  +  M).  tL  ♦  M  ♦  1 ) . K,  with  K  then  reset  to  K  ~  K  ~  M  as  noted  in  6*7.  The  re¬ 

placement  begins  at  L  rather  than  L  +  1  because  (A-l)  may  also  be  satisfied  by  (Xj ,  y^) and  some 
point  (xw .  yw ).  where  m  <  L.  Hence,  at  this  stage.  V  is  reduced  and  closed-up  (CU)  for  further 
processing.  Control  is  retunred  to  6-2  and  the  search  continues  through  the  updated  V-array, 
starting  with  i  a  k.  or  i  “  2  if  k  “  1  for  more  “duplicate"  points,  64. 

If  at  6-8.  IEN  -  K.  then  wc  must  have  k  =  1  (1ST  =  l)  when  ICV  >  0.  since  Vj  and  vK  are 
always  the  same.  Thus,  in  this  case.  6  proceeds  from  6-1 2  to  EXIT.  6-14.  and  return  of  control  to 
P-7.  If.  on  the  other  hand.  K'V  <  0.  i.e..  processing  of  V  is  from  N  -  I  to  l,  6*5,  and  IEN  -  K  at 
6-8,  and  1ST  =  j  /■  1  at  6-1 2.  then  an  element  has  been  isolated  which  is  specified  by  (j.  j  +  I , . . . .  K). 
Consequently  6-7  can  be  bypassed  with  V  reduced  and  CU  by  simply  resetting  xK  -*  xlST.  yK  -*■ 
Visr  ■  a,'d  1ST  -*  K  at  6-1 1 .  Figures  40.  52  contain  examples  of  where  this  would  occur. 

SORT  II.  7.  is  now  considered  in  more  detail  with  the  aid  of  its  flow  chart  7  (page  A- 19). 
Recall,  when  the  0-option  is  used,  tlwt  this  subroutine  is  called  by  SORT  1  to  delete  nodes,  from  an 


$Wo  use  p  ami  a  fur  an  isolated  dement  aiu!  retain  P  and  A  for  P(U)  and  A(U). 
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isolated  element  of  {S}  which  contains  SCP.  They  are  detected,  within  rounding  error,  by  testing 
the  inequality  (see  page  28), 


(A-2)  |si  -  Isin  Ad?j  <  to,  to  =  7(— 14), 

where  s  can  be  determined  algebraically  from  the  3  points  specifying  the  angular  region  (see  (47) 
on  page  28;  7-3,  12,  17).  Recall  -dso,  that  if  (A-2)  holds,  then  it  is  possible  |  AO  |  is  near  zero  rather 
than  n.  In  this  case,  although  the  angular  region  is  well  defined,  (WD),  the  vertex  node  is  deleted 
since  the  angular  region  does  not  contribute  to  p  or  a.  Angular  region  a,  in  Figure  32  is  an  example. 

It  is  assumed  now  that  an  isolated  element  S  of  {S}  is  available  through  the  decomposition 
ot  11  by  SORT  1.  We  assume  S  is  specified  by  an  array  r  of  M  coordinate  points.  Two  integer¬ 
valued  parameters  k  and  m  are  introduced  in  7-2  with  k  =  1 ,  m  =  2.  Parameter  k  is  associated  with 
the  vertex  point  (k)  of  the  angular  region  ak  under  consideration.  The  parameter  m  refers  to  the 
point  of  ak  following  ( k).  It  is  denoted  by  (m).  Initially  m  is  set  to  k  +  1,  7-2,  7-9.  If,  however, 
specified  by  (k  -  1,  k,  m)  subtends  an  angle  Af)  such  that  (A-2)  holds,  then  m  takes  successive 
values  above  k  +  1  until  (A-2)  is  not  satisfied  or  m  =  M,  74. 7-11, 7-16, 7-18. 

The  quantities  u,  v.  1)1  and  w,  z,  Dj  are  computed  initially  at  7-2  and  7-3.  respectively.  (The 
quantities  D,  and  !)>  are  defined  on  page  39.)  Then  (A-2)  is  checked  at  7-3  for  u( .  If  it  holds, 
then  m  -  m  +■  I,  74,  with  a  return  to  7-3  to  compute  new  values  of  w,  z,  Dj.  This  is  continued 
until  (A-2)  does  not  hold  or  m  ~  ,M.  If  tu  =  M  a  return  is  made,  7-5,  to  SORT  l.  box  13,  with 
NUMl  =  2.  Hence,  ptS)  =  0  for  this  particular  isolated  element  S.  since  it  is  a  straight  line  within 
the  tolerance  of  (A-2). 

Assuming  this  does  not  occur.  7  proceeds  to  7-7  with  ’  «  2  and  a  query.  Is  m  =  2?  If  in  2. 
then  points  tM ).  ( 1 ),  ....  (tn  ”  1)  were  found  to  be  colincar.  i.c. ,  each  3  successive  points  generate 
an  angular  region  for  which  (A-2)  holds.  In  this  case,  the  original  array  r  is  reduced  and  CU  by 

replacing  elements  of  r  starting  at  (t)  by  elements  (in  -  1).  Cm  -  2) . (M)  and  M  is  reset  to 

M  a  M  -  tm  ••  2 1.  7-8.  the  program  proceeds  to  7-9,  where  k  -  2.  m  a  3.  and  new  values  of  w,  t, 
and  Dj  arc  computed.  Then  7  would  proceed  to  7-t  2. 

tf  m  *  2  at  7-7,  then  is  WD  and  ?  proceeds  to  7-6,  without  disturbing  r.  with  m  -  3. 
Proceeding  to  7-12.  new  values  of  u,  v,  Dt  are  computed,  where  as  noted  above  k  -  2.  m  -  3. 
(Observe  that  at  this  stage  w,  z.  D;  from  7*3  are  based  on  k  s  1.  m  -  2.  and  therefore  have  the 
correct  values  for  looking  next  at  <u  )  lire  program  is  now  set  to  look  at  Uj.  where  the  subscript 
refers  to  clement  1 2)  of  the  updated  r  array. 

At  7-12  (A-2)  is  checked.  If  it  holds,  m  a  m  +  l  in  7-18  3m!  a  return  is  made  to  7-12.  where 
new  values  of  u,  v,  I);  arc  computed  and  (  A-2)  is  checked  again,  litis  is  continued  until  (A-2)  is 
not  satisfied  or  rn  =  M  ♦  1.  If  (A-2)  is  not  satisfied  for  some  m.  3  <  m  <  M.  then  from  7-12,  the 
program  proceeds  to  7-13.  If  in  a  k  +  1.  then  e,  is  WD.  no  alterations  are  made  lor.  k  sk  r  l. 
and  if  k  <  M.  the  program  goes  from  7-16  to  7-15,  where  Dt ,  u.  v  arc  used  for  the  new  values  of 
l)J.  vv.  i.  respectively,  m  -  k  +  l.  and  a  return  is  made  to  7-12.  lire  angular  region  a3  is  now 
investigated  as  was  done  previously  with  a  j.  If  k  —  M.  then  7  goes  from  7-16  to  7-17  io  process  aM. 
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If  at  7-13  m  #  k  +  1,  then  elements  of  r  starting  at  (k)  are  replaced  by  elements  (m  -  1), 
(m),  ....  (M),  with  M  reset  to  M  =  M-(m-k-  1),  7-10.  The  program  proceeds  to  7-11,  where 
k  =  k  +  1.  If  k  <  M,  7  returns  to  7-9  and  is  ready  to  look  at  the  next  angular  region.  If  k  =  M, 
then  7  proceeds  from  7-11  to  7-14  to  treat  the  last  angular  region  aM. 

An  answer  of  no  to  the  query  at  7-18  implies  ak  is  made  up  of  points  (k  -  1.  k.  M ; ;  conse¬ 
quently  all  the  points  (k  -  1),  (k),  ....  (M)  taken  as  successive  triplets  (k  -  1,  k,  k+  1),  (k-  l, 

k,  k  2) .  (k-1,  k,  M)  are  SCP,  i.e.,  satisfy  (A-2).  Therefore  points  (k  +  l ),  ....  (M-l)are 

ignored,  with  the  Mlh  point  replacing  the  kth  point  and  M  reset  to  k,  7-19.  It  remains  to  process 
aM  .  For  this,  we  go  co  7-14. 

Note  that  when  7  goes  from  7-16  to  7-17  to  compute  w.  z,  D,  foraM  that  u,  v,  are  already 
available  from  processing 

If  <iM  satisfies  (A*,!),  7-17  then  (M)  is  dropped  from  t  by  setting  M  =  M  -  l ,  7-20  and  control 
is  returned  to  SORT  I.  If  (A-2)  does  not  hold,  then  control  Is  returned  directly  to  SORT  I.  The 
final  coordinates,  as  contained  in  the  r  array,  at  exit,  and  the  number  of  them  are  specified  on  the 
flow  chart  of  7  as  output  x,  y,  M,  respectively. 

By  processing  polygons  11,  and  IK  of  Figures  32  and  33,  respectively,  a  more  detailed  descrip¬ 
tion  of  SORT  l  and  SORT  ti  is  given.  We  assume  the  (1-option  numbering  scheme,  in  order  to  bring 
SORT  11  into  play  for  I!,.  Also,  II,  and  IK  are  processed  in  the  order  of  increasing  numbered 
nodes,  starting  with  node  1.  Thus  1CV  =  3  (sec  P-7,  page  ( A-l 7)K  The  descriptions  are  presented 
in  tabulated  form  in  the  same  way  as  was  done  for  SORT  lit  (page  34 ).  Each  node  in  the  tabulation 
will  be  identified  by  its  number  in  the  original  V-array  specifying  the  given  polygon. 

The  first  column,  on  nage  A  12.  contains  the  value  of  N,  the  number  of  elements  fit  V,  when 
SORT  I,  6.  is  involved,  and  it  contains  the  value  of  M,  the  number  of  elements  in  the  r  array  when 
SORT  II.  7.  is  operating.  The  r  array  specifics  an  isolated  element  8  from  the  decomposition 
procedure  by  6.  The  second  and  third  columns  refer  to  integers  i  and  k.  and  k  and  m  of  the  preced¬ 
ing  discussions  on  6  and  7.  respectively.  The  fourth  column  displays  the  boxes  used,  by  their 
numbers  on  the  How  charts,  in  the  order  they  come  into  play.  C  olumn  four,  when  referring  to  6, 
also  shows  the  particular  S*  isolated  at  that  stage.  Column  five,  when  referring  to  7.  shows  the 
points  deleted  from  each  of  the-  8'  as  a  result  of  SCTTs).  The  numerical  data,  in  column  five,  pre¬ 
ceded  by  a  letter  is  associated  with  a  subsequent  column  headed  by  the  same  letter  which  shows  the 
reduced  CU  V  or  r-arrays  at  particular  stages  of  the  programs.  The  sixth  column,  headed  V, ,  refers 
to  the  original  V*array.  Tire  seventh  column,  headed  S, ,  refers  to  the  original  r  array  for  the  first 
Isolated  clement  S,.  Subsequent  elements  isolated  by  f>  have  their  initial  r  array*  listed  under 
columns  S:  and  S\  Columns  headed  (b).  (c).  etc.,  refer  to  the  reduced  compacted  arrays  as  deter¬ 
mined  by  7.  for  S'.  S2  and  S5.  For  example,  for  H,.  6  first  isolates  S1  given  by  (4.  6.  7.  8.  9,  lfi*. 

Then  S'  is  modified  by  deletion  of  (7.8).  lire  reduced  compacted  array  returned  to  VAl.R-7  is 
listed  in  column  beaded  (b).  Numerical  results  of  II,  arc  given,  following  its  tabulation,  at  the  end 
of  page  A- 13. 


For  Hj.  Figure  33.  SORT  1  decomposes  it  into  3  elements  of  the  class  >S? .  and  a  PAR.  speci¬ 
fied  by:  S'  *  (l.  2.  3.  4.  5.  (*.  7).  S2  «  (8.  9.  10.  1 1.  I2>.  S*  =  (12.  13.  14.  >5.  ’,6.  17.  !K>.  and 
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PROGRAMS  6  AND  7  FOR  11,  FROM  FLOW  CHARTS  (Continued) 
BASED  ON  FIGURE  32 


BOXES 

Compacted 

V  and  r-Arrays 

N 

i 

k 

SORT  I 

Points  Deleted 

S3 

(e) 

(0 

(g) 

(h) 

5 

4 

(V3):  1-3,  17-22 

1 

2 

2 

2 

2 

9 

10 

1 

2,4,6 

2 

3 

17 

17 

17 

S3  isolated  by  6 

3 

17 

18 

20 

20 

10,13 

17 

It 

19 

21 

21 

M 

k 

m 

SORT  II 

18 

19 

20 

22 

2 

9 

l 

2 

2,3 

19 

20 

21 

2 

9 

« 

3 

4,3 

20 

21 

22 

8 

2 

3 

7,  8,  9, 12 

(e);  1 

21 

22 

2 

8 

2 

4 

18, 12,  13 

22 

2 

7 

3 

4 

10, 11,9.  12 

(0:  3 

23 

7 

3 

5 

18,12 

7 

3 

6 

18,12,13 

5 

4 

5 

10,11.9,  12,  13 

(g):  18,19 

5 

5 

5 

16, 17 

4 

5 

5 

20,21 

(h):  22 

N 

i 

k 

SORT  1 

'  '1 

4 

5 

1 

9, 8,12,14 

Exrnop-7 

S4  Ml,  18,  19).  It  is  worth  noting  that  although  the  0-option  was  used,  SORT  11  is  not  needed. 
Tills  is  so  because  S*.  S3  and  S3  are  actually  in  (S)  and  S4  has  only  three  points  with  zero  area 
CNUM1  <  2  in  6-6).  Consequently,  SORT  11  can  be  bypassed  by  setting  ICV  =  2.  The  tabulation  for 
Il2  is  given  with  ICV  =  2. 

Hie  3  final  polygons  resulting  from  the  decomposition  of  II,  by  SORT  I  and  the  removal  of 
PAR(s)  by  SORT  II  are  listed  under  columns  (b),  (d)  and  (it).  VALR-7  yields  values  of  p  and  a 
for  cadi  of  these  polygons,  namely,  p(S‘)  =  -.7078  0769,  a(S')  =  -25,  p(S2)  =  .0268  8323, 
a(§J)  «  22,  p(S3)  ■  .0125  8574,  a(S3) «  13.5.  Tlius  P(I1,)  =  -.6683  3872,  A(I1,)  «  10.5.  It  is 
interesting  to  note  tliat  A  >  0  but  P  <  0.  In  Figure  53,  A  <  0  and  P  >  0. 
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PROGRAM  6  FOR  n2  FROM  FLOW  CHARTS 
BASED  ON  FIGURE  33 


BOXES 

Compacted  V-Arrays 
and  r-Arrays 

N 

i 

k 

SORTI 

Points  Deleted 

Vi 

b 

B 

S2 

B 

S3 

B 

18 

7 

1 

3,2,4 

1 

i 

i 

8 

l 

12 

l 

S1  isolated  by  6 

2 

2 

8 

9 

12 

13 

18 

(M  =  6 

6,10,9 

3 

3 

9 

13 

14 

19 

12 

2 

1 

8,7 

2-7  of  Vj 

4 

4 

11 

14 

15 

6 

2 

2,4 

5 

5 

11 

12 

15 

16 

(M  -  4 

S2  isolated  by  6 

6 

6 

12 

16 

17 

6,  10,9 

7 

7 

13 

17 

18 

8 

3 

1 

8,7 

8-11  of  V2 

8 

14 

18 

8 

2 

2,4 

9 

15 

19 

( 

M-6 

i 

S3  isolated  by  6 

10 

16 

6,  10.9 

11 

17 

2 

3 

1 

8,7 

12-1 7  of  Vj 

12 

18 

M  =  2 

^46 

13 

19 

8.  12,  14 

14 

EXIT  TO  P-7 

15 

16 

17 

18 

19 

The  decomposition  of  Il2  (Figure  33)  by  SORT  I  results  in  3  polygons,  S*.  S2,  $3.  Since  ai! 
the  angular  regions  of  these  polygons  arc  well  defined,  SORT  It  is  not  needed.  The  final  array  V, 
consists  of  a  singular  angular  region;  for  this  region  the  program  proceeds  from  6-6  directly  to  6-8 
setting  p  =  0  and  then  exiting.  VALR-7  yields  the  values  pfS1)  =  .8308  6076,  a(S')  =  18; 
p(S2)  =  .5378  8935,  a(S2)  =  6.5 ;  p(S3)  =  -.4271  2530,  a(S,)  =  “4.0.  Hence  P(Ha)  =  .9416  2481 . 
A(II-j)  =  20.5,  (see  page  48). 

Letting  S  denote  the  element  shown  in  Figure  36,  we  list,  on  page  A- 16,  P(ak)  for  each  angular 
region  ak,  k  =  1,  2 . N(=22).  The  computations  were  carried  out  with  IOP  =  3  (from  P-7),  Lc., 


with  9-decimal-digit  accuracy  for  P(ak),  for  each  k  =  1,  2,  N.  The  program  was  run  with 
ICV  =  0,  i.e.,  P-7  called  only  VALR-7  to  evaluate  P(S).  Subsequently,  the  program  was  also  run 
with  ICV  =  2  and  -2.  Recall  that  when  ICV  =  2,  P-7  calls  SORT  I  which  decomposes  n  (=S  here) 
into  a  set  of  simple  disjoint  polygons  {S1},  (in  this  case);  SORT  I,  in  turn,  calls  VALR-7  to  evaluate 
p(S‘)  for  each_ isolated  element,  S',  of  the  decomposition.  The  decomposition  by  SORT  I  starts  at 
point  (1)  of  S,  andjS  is  processed  from  (2)  to  (N),  sequentially.  When  ICV  =  -2,  the  procedure 
starts  at  point  N  of  S  and  carries  out  the  decomposition  in  the  backward  direction,  i.e.,  sequentially 
from  (N  -  1)  to  (1).  Observe  in  the  tabulation  that  the  decompositions  with  ICV  =  2,  and  ICV  =  -2, 
are  not  the  same,  although,  of  course,  the  final  results  for  P(S)  and  A(S)  must  be  identical  within 
the  accuracy  specified. 

The  first  column  of  the  tabulation  shows  the  node  number  of  S;  the  second  and  third  columns 
give  the  xy-coordinate  values  associated  with  the  node  number.  The  fourth  column,  headed 
ICV  =  0,  lists  P(ak)  for  each  node  number  (k)  of  the  first  column,  k  =  1,  2,  ....  N.  Summing  the 
P(ak),  and  using  (A-3),  below  gives  P(S)  beneath  columns  1-4.  The  next  two  columns  refer  to 
finding  p(S)  with  ICV  =  2.  The  fifth  column  contains  the  node  numbers  for  each  isolated 
S',  i  =  1,  ...,  6;  the  next  column,  headed  ICV  =  2,  contains  the  P(ak)  associated  with  the  node 
numbers  of  a  particular  S‘.  At  the  end  of  the  listing  for  each  S*,  p(S*)  is  given.  For  example, 
S3  is  specified  by  (13,  14,  15,  16);  the  value  of  P  for  the  angular  region  of  S3  at  node  13  is 
8.1418  04138  X  IQ"1.  The  '•  alue  of  p(S3)  =  7.3866  73215  X  10“2.  The  7<»  and  8“>  columns 
refer  to  nodes  and  corresponding  angular  regions  with  lf'V  -  -2.  For  example,  S4  in  this 
case,  is  specified  by  (4,5,6,7,12,13,17,18)  with  P(ti7 )  of  S4,  with  ICV  = -2,  given  by 
-2.0268  86540  X  lQ->  and  p(S4)  =  -9.1654  62410  X  I ()- 1 .  The  results  were  checked  by  using 
an  independent  decomposition  procedure,  with  Dreamer's  me* hod  (2),  described  on  page  47. 

It  was  shown  in  Section  III,  ((24)  and  (26))  that 

(A-3)  P(S)  =  l  -  £  P<ak>.  if  A(S>  >  0, 

(A*4)  P(S)  =  ~!  -  £  P(ak).  if  A(S)  <  0,  (See  pages  10-12) 

We  note  that  for  ICV  =  -2.  a(S' ).  a(S4)  are  negative,  (their  values  are  given  in  the  lower  right-hand 
corner  of  page  A*|6)>  i.e..  S1  and  S4  arc  negatively  oriented,  (NO),  and  therefore  p(S* )  and  p(S4) 
are  also  negative.  Note  again,  that  the  decompositions  for  ICV  =  2  and  ICV  =  -2  are  different. 
ICV  a  0  can  be  used,  because  §  has  no  SAR.  It  is  generally  preferred  when  no  SAR's  occur  for  S, 
since  it  does  not  use  SORT  l  nor  SORT  111  and  is  therefore  more  efficient. 

Summarizing  here,  we  can  say  that  (P-B)  is  significantly  better  than  (P-A)  for  complex 
SI  polygons  in  the  following  ways: 

(1)  Great  care  must  be  exercised  when  using  the  0-option  with  (P-A)  as  the  example  in 
Figures  A-5  and  A-6  shows.  Figure  58  is  another  example,  where  the  numbering  shown 
while  appropriate  for  (P-B),  since  VALR-2  handles  any  polygon,  is  inadequate  for 
(P  -  A)  for  the  same  reason  as  for  Figure  AT*. 

(2)  Program  B,  based  on  (P-B).  is  generally  more  efficient  tlian  the  Program  A  based  on 
(P-A).  i.e.,  VALR-7  with  SORT  I  and  11.  because  often  fewer  points  are  needed  to 
specify  11  as  in  Figure  58.  and  in  addition.  (P-A)  uses  an  NJ  process  to  decompose 
II  (SORT  1).  which  is  relatively  slow. 


TABULATION  OF  RESULTS  FOR  FIGURE  36 


1CV  =  0,  P(ak) 


1 

-5 

-5 

2.8665  15719  (-7) 

2 

5 

-5 

2.8665  15719  (-7) 

3 

5 

5 

2.8665  15719  (-7) 

4 

-5 

5 

9.9986  26366  (-1) 

5 

3 

3 

-1.3449  27576  (-4) 

6 

5 

-5 

-13736  33819  (-4) 

7 

-3 

-3 

-9.9864  95777  (-1) 

8 

3 

-3 

13480  75949  (-3) 

9 

3 

3 

9.9069  82439  (-1) 

10 

2 

-8.5731  82606  (-3) 

11 

-3 

-3 

-2.0268  94695  (-1) 

12 

-1 

0 

1.1979  54136  (-1) 

13 

2 

-1.8536  42054  H) 

14 

0 

1 

1.0265  18201  H) 

IS 

3 

3 

9.3010  33998  (-3) 

16 

-2 

-6.2292  $5125  (-4) 

17 

V.1 

3 

1.6245  62987  (-5) 

18 

-5 

5 

8. 2758  92935  (-2) 

19 

-3 

2 

8.9450  71S43(-1) 

20 

* 

2 

-9885?  13286  (-1) 

21 

-3 

-3 

-1.3654  79275  (-4) 

22 

-5 

S 

2.8665  157 19  (-7) 

»'(S)  -  1.6393  S343JM) 


tcv-o 

*S)-  I  -  Vl'fev).  A(S)«  $6.0 

I 


ICV  -  2 

t  i 

HS)  -  £  P<S‘).  A(S)  -2]  a(S‘) 
»  1 


tCV"  -2 

p(5)  -  2]  A(2)-V'a(st) 


ICV  =>  2,  P(ak) 


2.8665  15719  (-7) 
2.8665  15719  (-7) 
9.9986  26366  (-1) 
-1.3449  27576  (-4) 
2.7128  28364  (-4)* 


1.3496  06848  (-3) 
1.3480  75949  (-3) 
9.9069  82439  (-1; 
-8.5731  82606  (-3) 
1.5177  25594  (-2) 


8.1418  04138  (-1) 
1.0265  18201  (-1) 
9.3010  33998  (-3) 
7.3866  73213  (-2) 


2.2232  56327  (-2) 
1.6245  62987  (-5) 
8.2758  92935  (-2) 
8.94 SO  71843  (-1) 
4.8S07  74211  (-4) 


7.9730  92908  H) 
1.1979  54136  (-1) 
9.0285  63470  (-3) 
7.3866  73215  (-2) 


231665  15719  (-7) 
9.9986  16366  (“1) 
-I  J449  27576  H) 
2.6665  157191-7) 
2.71 28  28364  (-4) 


ICV  =  -2,  P(ak) 


-91724  10707  (-1) 
8.9450  71843  (-1) 
-9.6857  13286  (-1) 
-1.3654  79275  (-4) 
-8.5582  37140  (-3)* 


8.1418  04138  (-1) 
1.0265  18201  (-1) 
9.3010  33998  (-3) 
7.3866  73215  (-2) 


1.3496  06848  (-3) 
1.3480  75949  (-3) 
9.9069  82439  (-1) 
-8.S731  82606  (-3) 
1.5177  25594  (-2) 


-1.3736  33819  (-4) 
-1.3449  27576  (-4) 
-1.3736  33819  (-4) 
-2.0268  86540  (-1) 
1.1979  54136  (-1) 
-1.6754  46491  (-4) 
1.6245  62987  (-5) 
“9.1654  62410  (  -1) 


2.866$  15719  (-7) 
28665  15719  (-7) 
28665  15719  (-7) 
2.8665  15719  (-7) 
9.9999  ,88534  (-1) 


*(S‘)  «  20.  *(8J)  *  6.  a(S*)  -  38  |  a(S')  -  -7.S.  a(SJ)  «  3.5,  a(Sj)  ■  6 
a(S«)  -  3.  a(S*)  -  38.  a(S*)  -  20  I  a(S4)  -  -46.  a(S>)  •  100. 


*No4c:  Set  (A-3)  an*!  (A-4)  at  page  (A-l  S), 


'When  SORT  I  is  called,  il  processes  the  arbitrary  polygon  II.  decomposing  it  into 
elements  S*.  . ...  S}  with  aoption,  or  S*  ....  Sf  with  jjopiion,  VALR-7  is  called  to 
compute  P(S* )  or  PCS')  and  sums  the  results  to  get  P(fl).  In  case  of  S’.  SORT  1  calls 
SORT  II  to  remove  SCP.  Details  in  Flow  Chart  7. 
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EVERY  SIMPLE  POLYGON  CONTAINS  AN  INTERIOR  DIAGONAL 
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APPENDIX  B 

EVERY  SIMPLE  POLYGON  CONTAINS  AN  INTERIOR  DIAGONAL 


By  an  interior  diagonal,  we  mean  the  open  segment  (L)  of  the  closed  line  segment  [L]  extend¬ 
ing  from  some  vertex  (k)  to  some  nonadjacent  vertex  (j)  of  the  simple  polygon  S,  such  that  (L)  is 
entirely  in  the  interior  of  S. 

The  proof  given  in  (4,  p.  17,  Lemma  2]  is  not  correct.  Knopp’s  proof  is  repeated  here  with  a 
counter-example.  An  argument  to  correct  the  proof  is  then  given. 

Let  a  straight  line  which  does  not  intersect  or  meet  S  be  translated  parallel  to  itself  toward  the 
polygon  until  they  meet.  Then  the  line  necessarily  contains  a  vertex  A  of  S  with  the  interior  angle 
of  A  less  than  two  right  angles.  Let  B  and  C  denote  the  adjacent  vertices  to  A.  Then  one  of  the 
following  is  true: 

( 1 )  BC  is  a  diagonal  lying  in  the  interior  of  S. 

(2)  There  is  at  least  one  vertex  of  S  on  the  (open)  segment  BC  (let  one  of  these  vertices  be 
denoted  by  V)  but  no  vertex  in  the  interior  of  triangle  ABC,  (A  ABC). 

(3)  There  is  at  least  one  vertex  of  S  in  the  interior  of  A  ABC- 

If  (1)  is  true,  there  is  nothing  further  to  show.  If  (2)  holds,  then  AV  is  an  interior  diagonal  of  S. 
If  (3)  is  true,  let  a  point  X  move  from  B  to  C  along  BC  until  AX  encounters  a  vertex  or  vertices  of  S 
in  the  interior  of  A  ABC.  If  V  denotes  tliat  one  of  these  vertices  wluch  is  neatest  to  A,  then  AV  is  a 
diagonal  interior  to  A. 

The  proof  of  part  (3)  is  not  correct.  This  is  easily  seen  from  the  figure  below.  The  vertex 
nearest  to  A  in  A  ABC.  following  Knopp.  is  D,  but  the  line  AD  contains  points  outside  S.  The 
proper  vertex  to  have  chosen  was  li  which  is  in  A  ABC.  but  farther  from  A  than  D. 


A 


B-3 


The  proof  is  easily  corrected  as  follows:  Starting  at  A,  move  an  open  segment,  which  extends 
from  AB  to  AC,  parallel  to  BC  and  towards  BC  until  one  or  more  vertices  of  S  are  met.  If  there  is 
more  than  one  such  vertex,  choose  any  one  and  call  it  V.  Vertex  V  has  the  property  that  no  other 
vertex  of  S  in  A  ABC  has,  a  greater  (perpendicular)  distance  from  BC.  Now  suppose  AV  is  not  an 
interior  diagonal  of  S.  Then  there  exists  a  point  (z)  where  AV  intersects  another  side  of  S,  say  side 
(k,  k  +  1 ).  Point  (z)  cannot  be  a  vertex  by  the  way  V  was  chosen.  Now  either  vertex  (k)  or  (k  +  1 ) 
must  have  at  least  as  great  a  distance  from  BC  as  (z).  Say  it  is  (k).  But,  since  S  is  simple  (k)  must  be 
in  the  interior  of  A  ABC.  This  contradicts  the  way  V  was  chosen.  Hence  AV  must  be  an  interior 
diagonal  of  S. 

This  result  is  used  on  page  1 1  and  in  Appendix  D. 
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APPENDIX  C 

AN  ALTERNATIVE  METHOD  TO  FIND  P  FOR  SIMPLE  POLYGONS 


1 


R’ 


fj 

$ 

V| 

C 

f 
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At  an  early  stage  of  our  studies,  we  developed  a  method  to  compute  the  P-function  over  a 
simple  polygon  by  using  a  program  already  available,  which  computed  P  for  convex  polygons,  [2]. 
To  put  it  another  way,  an  automatic  procedure  was  set  up  to  represent  any  simple  polygon  by  a 
finite  set  of  convex  polygons.  Realization  of  our  working  program  required,  in  addition  to  a 
program  for  computing  P  for  convex  polygons,  a  program  to  determine  the  convex  hull  of  a  finite 
point  set  in  the  plane.  Such  a  program  was  available  from  previous  work,  [  1  ] .  By  the  convex  hull 
C(ZN)  of  the  point  set  ZK  =  {(xj,  yj),j  =  1,  2, N},  we  mean  the  smallest  convex  polygon  which 
contains  all  of  ZN.  The  vertices  of  C(ZN)  are  in  ZN. 

A  simple  polygon  S  is  shown  in  Figure  C-l.  We  set  forth  the  procedure  by  applying  it  to  this 
polygon.  It  will  be  apparent  to  the  reader  that  any  N-sided  simple  polygon  can  be  handled  in 
the  same  way. 

Procedure: 

(A)  Find  the  convex  hull  C  of  S.  We  obtain 

C  =  (1,2,  3,  6,  7,  13,14). 

The  P-function,  P(C),  for  C  is  computed  by  the  program  for  evaluating  P  for  convex 
polygons.  Clearly  since  S  and  C  are  positively  oriented,  PO,  we  have 

(C-l)  0  <  P(S)  <  P(C). 

(B)  The  set  of  vertices  of  C  is  searched  to  determine  which  vertices  of  S  are  missing  between 
adjacent  vertices  of  C.  Obviously,  vertices  (4)  and  (5)  between  (3)  and  (6),  and  vertices 
(8),  (9),  (10),  (11),  (12)  between  (7)  and  (13)  are  missing  from  C.  In  this  way,  we 
isolate  2  simple  polygons  from  C,  namely 

Sj  -  (j,  4,  5,  6,  3) , 


Figure  C-l.  A  1 3-Sided  Simple 
Polygon,  S 


S2  =  (7,8,9,10,11,12,13,7). 
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(C)  The  convex  hull  of  S5  is  found  to  be  identical  to  S3.  (When  this  occurs,  that  convex 
hull  requires  no  further  processing.)  We  note  C;  is  negatively  oriented  so  that  P(Cj)  <  0. 
The  convex  hull  C2  of  S2  is  then  determined  to  be 

C2  =  (7.9,  10,  13,7), 

as  indicated  in  Figure  C-2  by  the  dotted  lines  and  the  line  segment  (9,  10).  C2  is  also 
NO,  and  hence  P(C2)  <  0.  Therefore,  consideration  of  C,  Ct,  C2  shows 

(C-2)  P(S)  >  P(C)  +  Pfq)  +  P(C2). 

(D)  Two  simple  polygons,  both  PO,  are  obtained  from  C2,  by  noting  the  missing  vertices,  as 
explained  in  (B),  namely 

S3  =  (7,  8,  9,  7),  S4  =  (10,11,12,13,10). 

(E)  The  convex  hull  C3  for  S3  is  again  S3.  Thus,  C3  requires  no  further  processing.  Since 
C3  is  PO,  we  obtain  P(C3)  >  0.  Next,  the  convex  hull  C4  for  S4  is  found  to  be 

C4  =  (10,  12,  13,  10) 

'ch  is  also  PO.  Thus  P(C4)  >  0  and  we  obtain,  adding  P(C3)  and  P(C4)  totheright- 
nand  side  of  (C-2), 

(C-3)  P(S)  <  P(C)  +  P(C,)  +  P(C2)  +  P(C3)  +  P(C4). 

(F)  Finally  wc  isolate  Ss  from  C4> 


Ss  =  (10, 11,12, 10). 

Its  convex  hull  Cs  is  identical  to  it.  Observing  that  Cs  is  NO  and  that  P(CS)  <  0,  we 
obtain  the  final  results  by  adding  P(CS)  to  the  right-hand  side  of  (C-3),  namely,  with 
C0  sC, 


0  c.  ■  s. 

t-0 


« 
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Figure  C-2.  Convex  Hull,  C2,  of 
Sj  “(7,8, 9, 10, 11, 12, 13, 7) 


o 


C4 


Although  no  proofs  are  given  the  procedure  can  be  put  on  a  rigorous  basis  by  induction 
type  arguments. 

This  method,  cal!  it  (0)  (for  old)  is  much  slower  in  general  than  the  procedure  described  in 
Section  III,  call  it  (N)  (for  new).  This  is  so,  because  (in  addition  to  finding  convex  hulls)  the  time 
consuming  computation  is  finding  P  for  an  angular  region  by  VALR-2  or  VALR-7.  By  our  present 
procedure,  (N),  we  require  the  evaluation  of  P,  in  the  example  of  Figure  C-l,  for  13  angular  regions, 
whereas  by  the  method  of  this  appendix,  (0),  we  require  6  angular  regions  for  P(C),  4  foi  P(Cj ), 
4  for  P(C2),  3  for  P(C3),  3  for  P(C4)  and  3  for  P(CS)  for  a  total  of  23  angular  regions. 

The  method  (0)  does  have  an  advantage  in  that  P(S)  is  alternately  bounded  above  and  below 
with  improved  bounds  on  each  cycle  of  positive  and  negative  contributions  to  estimating  P(S). 
By  a  cycle,  we  mean  a  stage  in  (0)  where  each  convex  hull  obtained  is  of  the  same  orientation. 
The  first  cycle  occurs  with  C  is  found.  It  is  PO.  At  the  second  stage  Ct  and  C2  are  found  and 
both  are  NO.  The  third  stage  is  manifested  by  the  appearance  of  C3  and  C4,  both  PO.  The 
fourth  and  final  stage  is  when  Cj  is  found.  It  is  NO.  The  improved  bounds  may  allow  the  calcula¬ 
tion  for  P($)  to  be  terminated  early.  Indeed,  if  at  the  end  of  any  cycle  of  NO  convex  polygons,  the 
last  denoted  by  Cs ,  the  quantity 


J 

(C-5)  £  P(Cj) 

i*0 

is  greater  than  1  -  e,  then  P(S)  »  1  within  e;  if  at  the  end  of  any  cycle  of  PO  convex  polygons,  the 
last  denoted  by  Cj ,  the  quantity  corresponding  to  (C-5)  is  less  titan  s,  then  P(S)  -  0  within  e. 
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We  first  show  an  expression  for  the  area  of  a  simple  polygon  in  terms  of  vectors.  There  is 
nothing  new  about  the  result,  but  it  is  not  as  easily  available  as  one  would  expect,  [6], 
Subsequently,  by  some  straightforward  algebraic  manipulations  we  obtain  an  expression  for  the 
area,  A,  which  leads  to  a  very  efficient  machine  computation  for  the  area  of  a  simple  polygon,  A(IT). 

Let  S  denote  a  simple  polygon  with  its  vertices  numbered  in  the  natural  order  from  1  to  N, 
such  that  in  tracing  S  continuously,  the  interior  of  S  is  always  on  the  left.  We  say  S  is  positively 
oriented  (PO)  in  this  case.  The  classical  vector  expression  A,  for  which  I  A(S)|  =  A(S),  is  then  given 
by  a  sum  of  vector  cross-products. 

1  N 

(D-D  A(S)  =- £  (z,-Z)X(zl+1 -Z).  2j\j  + 1  =  Zj , 

i*  1 

where  5,-2  denotes  the  vector  from  2  to  5,,  with  Z  fixed,  but  arbitrary,  and  A(S)  >  0  using  the 
right-hand  rule  for  cross  products. 

In  order  to  establish  (D-l)  for  simple  polygons,  we  use  induction,  and  the  result  of  Appendix  B 
that  for  any  simple  polygon  S  there  exists  a  diagonal  between  two  vertices  of  S  which  is  entirely 
in  S.  We  also  need  the  fact  that 

(D-2)  (z,  - 2)  X  (Zj  -  Z)  =  -(Zj  -  2)  X  (5,  -  2). 

Certainly  for  N  =  3.  S  a  triangle,  (D-l)  holds.  Now  assume  (D-l)  holds  for  all  simple  1*0 
polygons  with  no  more  than  N  -  l  vertices.  Let  S  denote  a  simple  polygon  of  N  vertices,  PO.  By 
the  result  of  Appendix  B,  them  exists  a  diagonal  from  vertex  j  to  vertex  j  +  k,  k  >  1 ,  which  is 
entirely  within  S,  except  for  its  end  points  at  vertices  j  and  j  +  k.  This  diagonal  divides  S  into  two 
disjoint  simple  PO  polygons,  except  for  the  common  side,  each  with  no  more  than  N  -  1  vertices. 
Hence,  by  the  induction  hypothesis,  (D-l)  holds  for  each  of  these  polygons,  call  them  S,  and  S2, 
where 


Sj  =  (1.  3, ....  j  -  1.  j.  j  +  k,  j  +  k  +  I. N,  I), 
S2  3  (j.  j  +  1 . j  +  k-  1,  j  +  k,  j). 

Hence 


(D-3) 

A(Sj ) 


l-i 


£  (z,-Z)X(5,*,  -Z)  +  (zrZ)X(zjHi-Z)  +  £  (5,-2)  X  (5,*,  -2) 


U«i 


i*j*k 


(D4) 


A(Sj)  =  ± 


1  ♦k-t 


y .  (^t  ~  zj  x  (z,», 


t*i 


-Z)  +  (*,♦„  -Z)X(zr2) 
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Since  Sj  and  S2  are  disjoint  and  PO,  we  have  by  adding  (D-3)  and  (D-4),  and  using  (D-2),  the 
expression  (D-l)  for  N-sided  simple  polygons. 

Above,  we  have  assumed  S  was  PO.  If  S  is  NO,  then  each  cross  product  in  (D-l)  is  reversed, 
and  by  (D-2),  A  is  given  by  (D-l)  with  a  minus  sign  attached. 

Now,  since  A(S)  is  a  continuous  function  of  the  coordinates  of  the  vertices  of  S,  (D-l)  also 
yields  A  for  polygonal  elements  of  {S},  such  as  in  Figures  35, 36, 45. 

Actually  (D-l)  holds  for  arbitrary  polygons,  i.e.,  elements  of  {II}.  This  follows  by  using  the 
above  results  with  a  theorem  given  in  [4,  page  15],  which  states  that  every  polygon  can  be 
decomposed  into  a  finite  set  of  polygons  in  (S). 

An  efficient  expression  for  computing  A  can  be  obtained  from  (D-l).  Since  Z  is  arbitrary, 
choose  Z  =  6.  Then  (D-l)  reduces  to 


(D-5) 


A(II)  =  Jj-  (ZjXzj+i). 
1 


In  component  form  we  have 


Z(XZi  +  ,  -*•  x,y,*, 


so  that  (D-5)  can  be  written  as 

1  N-> 

(D-6)  A(I1)  =  ^  ]T  (xtyi  +  ,  -xmy,)«  (xN*,,  yN  +  J)  s  (x,.  y,). 

“  i-i 

The  number  of  multiplications  can  be  halved  by  some  algebra.  From  (D-6)  take  the  second  product 
of  the  (i  -  l  )a  term,  x(ylw|,  and  combine  it  with  the  First  product  of  the  ilU  term,  x(y1+ ,  to  obtain 
xi(>Vi  "yi-j).  This  can  be  done  successively  for  each  i  °>  2,  3,  ....  N.  The  remaining  elements, 
namely,  the  first  product  of  the  first  term,  x,  y2  and  the  second  product  of  the  N,b  term,  x,yN  are 
combined  to  obtain  Xj  (y2  -  yN).  Thus  (D-6)  becomes 

i  N 

(D-7)  A(I1)  £  Xj(ym -y,_,).  yN*,  =  y,,  y0syN-  (See  (22)  and  (46)). 

i-i 

This  expression  appears  in  the  text  as  (22)  and  (46). 

A  Fortran  IV  listing  of  the  short  program  for  computing  A,  SMP-7,  is  given  in  Appendix  F, 
page  (F-37). 

The  area  of  II  in  the  wz-piane,  A(w,  z),  sec  page  i ,  is  given  by 
(D4J)  A(w,  z)  «  owot(  1  -pJ)‘/2A(m. 
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In  this  appendix,  we  list  the  constants  that  appear  in  the  programs  SORT  I,  SORT  II,  SORT  III, 
VALR-2,  and  VALR-7.  There  are  two  constants,  a  and  cj,  which  depend  only  on  the  characteristics 
of  the  computer  used.  They  are  set  at  5  X  10-d  and  7  X  10~d,  where  d  denotes  the  maximum 
number  of  decimal  digits  the  machine  uses  to  represent  a  real  number.  For  our  machine  d  =  1 4. 

For  VALR-2  and  VALR-7,  the  additional  parameters  that  appear  are  listed  for  4  levels  of 
accuracy,  i.e.,  3,  6,  9  and  12-decimal  digits.  The  last  is,  at  present,  not  incorporated  into  our 
programs,  but  it  would  be  easy  to  do  so.  The  values  for  all  the  parameters  follow. 

o  =  5  (- 1 4)  =  5  X  1 0- 14 ,  in  SORT  1 
=  7 (- 1 4)  in  SORT  11  and  SORT  III 

w  =  7 (—  1 4)  in  VALR-2  and  VALR-7 


ADDITIONAL  PARAMETERS  FOR  VALR-2.  OR  VALR-7 


Acc. 

e 

n 

“2 

(R/x/2)2  =  c2 

® 

2.54  (-4) 

2.02  (-7) 

1.22  (-2) 

5.625  (-5) 

6.962  (-2) 

6.05160 

2.57  (-7) 

2.08 (-13) 

1.23  (-4) 

5.700  (-8) 

6.990  (-3) 

12.60605 

© 

2.94  (-10) 

2.71  (-19) 

1 .34  (-6) 

6.512  (-11) 

7.311  (-4) 

19.201924 

1.00  (-13) 

3.17 (-26) 

6.58  (-9) 

2.225  (-14) 

5.111  (-5) 

26.103925 

e  “  6/>/i 

See  page  6 

1  t»j  =  y/v  c/2  See  (2 

.page  15) 

a,  -vc2  See  page  7  a4  =  otj j  Sec  page  29,  liq.  ( 1 2)  also, 

otj  =  (9a,),/J  Sec  page  7  K1 12  See  pages  6.  28  and  (2,  page  81 


The  first  column  of  the  table  labeled  Acc.  (for  accuracy)  lists  ®  (§)  referring  to 
3,  6,  9.  1 2-dccimal-digit  accuracy,  respectively,  for  the  probability  over  an  angular  region.  Pages 
arc  indicated  above  where  the  parameters  are  discussed  in  the  report. 

Tire  ininimax  coefficients,  ak.  for  approximating  erfc  (x)/z(x)  on  (0,  c(5)l  (see  page  6)  are 
given  below  for  the  four  accuracy  levels  associated  with  ®,  (B),  ©>  ©  as  noted  above.  They  were 
computed  by  a  double  precision  minimax  subroutine  utilizing  values  of  erfc  (x)  correct  to  18 
significant  digits  on  ( 1/2,  c(8))  and  values  of  erf  (x)  accurate  to  25  digits  on  |0, 1/2). 
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For  (A)  (Average  time  per  angular  region  =  7.8  X  10-4  sec) 


a0  =  .8857775 18572895D  +  00 
a2  =  .759305502082485D  +  00 
a4  =  .695232092435207D  -  01 


t  =  -.981 151 952778050D  +  00 
3  =  -.353644980686977D  +  00 


For  (g)  (Average  time 


per  angular  region 


1-1  X  lO-3  sec) 


a0  =  .8862264700 16632D  +  00 

a2  =  -885348820003892D  +  00 

a4  a  -4218211 97 160099D  +  00 

a6  =  .905057384 150449D  -  01 

a8  «  .430895 1689841 38 D  ~  02 


‘i  =  -.99995071 456 1036D  +  00 

,3  =  -.66061 1239043357D  +  00 

5  =  --222898055667208D  +  00 

7  =  -.2549061 11884287D  -  01 

9  =  --323377239693247D  -  03 


For  (£)  (Average  time  per  angular  region  =  1.3  X  IQ- 3  sec) 


a0  =  .88622692493 1465D  +  00 
a2  °  .886223733 186722D  +  00 
«  .442851899328568D  +  00 
a5  =  .1450600434030121)  +  00 
a8  a  .309 1 9929552 1 2 1 0D  »  01 
aio°  .324944543171 185D  -  O'* 
ai  2=  •  1 05787574480633D  -03 
a14“  .40833551 7232 165D  -06 


a,  =  -.99999989977625  2D  +  00 
a3  =  .6666266705109071)  +  00 

a5  =  -.2656382063660251)  +  00 
u7  =  -.7149098377998891)  -  01 
a9  *  -•!  1 2323532 1 4844 1 D  -  01 
an  =  ", 7042602433090961)  -  03 
at3  »  -.97 1864864 16046 ID  -  05 


For  (0)  (Average  time  per  angular  region 


ao  *  .8862 2692545  2593 D  +  00 
a2  «  .8862269227S6746D  +  00 
a4  •  .4431128680489191)  +  00 
a6  -  .1476871363219381)  +  00 
h  =  .3680328493508600  -  01 
aio«  .7 1 0292625 73405 2 D  02 
ai  2“  981I12629090333D  -  03 
an=  .789960968802448D  -  04 
ai6=  .283646635409.3221)  -  05 
al8»  .317679497040006!)  -  07 
a20=  4525343473373050  -  10 


1.5  X  10-3  sec) 


al  s  .9999999999485970  +  00 

aj  ■  6666666118666610  +  00 

a5  =  -  26666272909141 1 0  +  00 

a7  a  -.7613658558502920  01 

a9  *  -1671950968881830  -  01 

ali“  -.2781709329062240  -  02 
an*  -.3025886407521080  -  03 
al5=  - .  1 68685 1 8 1 7670460  -  04 
-.3583144669082900  0f> 

al9“  “  .  1 7544065 1 9404300  -  08 


Aveng,  time  per  angular  region  e  relcm  lo  die  avenge  computing  time  on  the  CDM  700  to  obtain 


$  r'.'<A^vyV..Vl»  X1/'  »r  - 


APPENDIX  F 

PROGRAM  LISTINGS  IN  FORTRAN  IV 
P-2,  VALR-2,  SORT  III,  VALR-7,  P-7,  SORT  I,  SORT  II,  SMP-7 
(Flow  cJiaris  on  pages  40-45  and  A-I7  to  A-I9) 


MASTER  SUBROUTINE  P-2 
(FLOWCHART  1,  page  40) 

P-2  is  used  for  computing  P(II)  over  an  Arbitrary  Polygon  II* 

CALL:  P-2  (x,  y,  N,  P,  ICV,  IND,  IOP,  A,  W), 
where: 

x  is  the  array  of  abscissas  of  the  numbered  points  of  n.  x  is  dimensioned  at  N  +  1 . 

y  is  the  array  of  ordinates  of  the  numbered  points  of  II.  y  is  dimensioned  at  N  +  1 . 

N  is  the  number  of  points  specifying  n,  except  if  N  =  1  when  the  IBND  over  an  angular 

region  is  computed.  Three  input  points  are  needed  when  N  =  1,  given  in  counterclock¬ 
wise  order,  with  the  vertex  at  point  one,  (see  pages  25,  27). 

P  is  the  location  where  the  value  of  P(IT)  is  returned. 

ICV  must  be  set  as  an  integer  by  the  user  according  to  the  list  below: 

ICV  =  0,  n  is  simple,  or  of  S  type  with  no  SAR(s)  (see  pages  12, 31).  VALR-7  used  alone. 
ICV  >  0,  II  is  arbitrary.  VALR-2  used  alone. 

ICV  =  -2,  II  is  of  S  type  with  possible  SAR(s). 

ICV  <  0,  i£-2,  II  is  arbitrary  with  PAR(s). 

IND  is  an  error  indicator.  Normally,  it  is  set  to  zero.  If  IND  =  2,  then  PAR(s)  have  been 
detected  by  either  VALR-2  or  VALR-7.  For  VALR-2,  (ICV  >  0,  ICV  <  0,  *  -2)  the  re¬ 
sult  for  POl)  is  acceptable.  For  VALR-7  (ICV  =  0,  -2)  however,  this  result  of  IND  =  2, 
means  the  value  for  P(II)  is  most  likely  wrong,  unless  N  «  1  VALR-7  is  not  to  be  used 
alone  where  SARfs)  are  a  possibility,  unless  N  =  I,  If  IND  =  3,  then  N  has  not  been 
specified  as  an  intogor  equal  to  one  or  greator  than  two.  Such  values  of  N  are  not 
allowed. 

IOP  Is  an  accuracy  parameter.  It  is  set  by  the  user  to  1,  2,  3  for  approximately  3,  6,  or 
9  decimal  digits  of  accuracy  in  P(ri). 

A  is  the  location  where  A(  II)  is  returned.  I A  i  gives  area  of  II,  (see  pages  9, 26). 

W  is  the  location  where  the  winding  number  of  II  is  returned.  It  is  computed  in  VALR-2 

and  takes  integer  values  (see  pages  18,  19).  W  is  defined  as  an  integer  variable.  It  is 
initialized  to  one,  and  is  only  computed  if  ICV  >  0  or  ICV  <  0,  =*  -2. 


*Sce  footnote  ! ,  page  I ,  for  definition  of  an  arbitrary  polygon. 
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SUBROUTINE  P2  (  X ,Y,NB,P, ICV, IND, IOP, A,KO  ) 
DIMENSION  X(l) / Y (1) 

IF  (  NB.NE. 2. AND.NB.GE. 1  )  GO  TO  3031 

IND=3 

RETURN 

3031  CONTINUE 
N=NB 
KO=l 

IF  (  ICV.EQ.O.OR.NB.EQ.l  )  GO  TO  3091 
IF  (  ICV.GT.O  )  GO  TO  3071 
CALL  SORT3  (  X,Y,N  ) 

IF  (  N.GT.2  )  GO  TO  3061 
P=0. 

A=0. 

IND-0 

RETURN 

3061  CONTINUE 

IF  (  ICV.EQ.-2  )  GO  TO  3091 
3071  CONTINUE 

CALL  VALR2  (  X , Y,N,P, IOP, A, IND,KO  ) 

RETURN 

3091  CONTINUE 

CALL  VALR7  (  X,Y,N,P, IOP, A, IND  ) 

RETURN 

END 
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SUBROUTINE  VALR-2 
(FLOW  CHART  2,  page  41) 

VALR-2  is  used  to  compute  P(I1)  when  II  is  arbitrary 


CALL:  VALR-2  (x,  y,  N,  P,  IOP,  A,  IND,  W), 
where: 

x  is  the  array  of  abscissas  of  the  numbered  points  of  II.  x  is  dimensioned  at  N  +  1 , 

y  is  the  array  of  ordinates  of  the  numbered  points  of  II.  y  is  dimensioned  at  N  +  1 . 

N  is  the  number  of  points  specifying  II,  except  if  N  =  1  when  the  IBND  over  an  angular 

region  is  computed.  Three  input  points  are  needed,  when  N  =  1,  given  in  counterclock¬ 
wise  order,  with  the  vertex  at  point  one,  (see  pages  25,  27). 

P.  A  are  the  locations  where  the  values  of  P(I1)  and  A(II)  are  returned, 

IOP  is  an  accuracy  parameter.  It  is  set  by  the  user  to  1,  2,  or  3  for  approximately  3, 6,  or 
9-decimal  digits  of  accuracy  in  P(II). 

IND  is  an  error  indicator.  Normally,  it  is  set  to  zero.  If  IND  =  2,  it  informs  the  user  that  II 

contains  a  PAR.  The  value  for  P(II)  is  acceptable,  if  IND  =  3,  then  N  has  not  been 
specified  as  an  integer  equal  to  one  or  greater  than  two.  Such  -alues  of  N  are  not 
allowed. 

W  is  the  location  where  the  value  of  the  winding  number  W  for  II  is  returned.  W  is  an 
integer  variable. 

This  routine  requires  computation  of  erf  (x)  and  erfc  (x)  which  are  defined  on  pages  S,  28  and 
29.  We  have 


ERI  .  (x)  =  erf(x),  ERFC  1  (0,x)  «  ertc(x), 

where  the  subroutine  listings  for  these  functions  are  given  on  pages  F-1 2  to  F-l  5.  They  are  identical 
to  the  NSWC(DL)  math  library  functions  ERF  and  ERFC  as  of  June  1980  which  are  based  on  the 
reference  below. 

Cody,  W.  J„  Rational  Chehyslm ■  Approxitmtlonr  for  the  Error  Function,  Mathematics  of 
Computation,  v.  23  (1969),  pp.  631-637. 
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SUBROUTINE  VALR2  {  X ,Y ,N, P, IOP, A, IND,KO  ) 
DIMENSION  X(l) ,Y (i)  ,G(2) ,H(2) ,RSQ(4) 
DIMENSION  E<5) ,E2(10) ,E3(15) 

DIMENSION  APH1{3) ,APH2(3) ,CST(3) 

DIMENSION  APH4 (3) ,  A3D8 (3) 


REAL  L 
REAL  ROM 

DATA  PI/3.1415  92653  5898  / 
DATA  TWOPI/6.2831  85307  17958 
DATA  ALNPI/1. 1447  29885  84940 
DATA  Cl/. 28209  47917  73877  / 
DATA  C2/. 15915  49430  91895  / 
DATA  TAU/7.E-14  / 

DATA  TAUSQ/4 . 9E-27  / 

DATA  (  E(I),I=1,  5)  / 

. 885777518572895E+0Q  , 
. 759305502082485E+00  , 
. 695232092435207E-01  / 
DATA  (E2(I),I=1,  10)  / 

. 886226470016632E+00  , 
. 885348820003892E+00  , 
. 421821197160099E+00  , 
. 905057384150449E-01  f 
. 430895168984138E-02  , 
DATA  (E3{I),I«1,  15)  / 

.88622692493 146 5E+0Q  , 
. 8Q6223733186722B+00  , 
, 442851899328569E+00  , 
. 14 506004 340 3014E+00  , 
. 309199295521210E-01  r 
.32494454317U85E-02  , 
. 105787574480633B-03  , 
.408335517232165E-06  / 
DATA  {  APK1 (I) , l»l,  3  )  / 
2.02E-7,2»08£-13,2.71E-19  / 
DATA  (  APH2(I) ,1*1, 3  )  / 
1.228-2,1.236-4,1. 34E-6  / 


/ 

/ 


-.9811519527780502+00  , 
353644980686977E+00  , 


- . 999950714561036E+00  , 
- .660611 23904 33 57E+00  , 
- . 222898055667208E+00  , 
- . 2549061118B4287E-01  , 
- .323 37723 9693247E-03  / 

- . 999999899776252E+00  , 
- .6666266 7 051090 7E+00  , 
- . 265638206366025E+00  , 
- . 7149Q9837799889E-01  , 
- . 112323532148441E-01  , 
- . 704260243309096E-03  , 
-.97 186486416046 IE-05  , 


DATA  (  APH4 (I) ,1*1,3  )  / 

.6962E-1,  . 6990E-2,  .7311E-3  / 

DATA  RTPII/. 56418  95835  4776  / 

DATA  {  RSQ  \'I) ,  1-1,3  )  / 
6.0516,12.60605  ,19.201924  / 

DATA  (  A3D8 (X) ,1=1,3  )  / 

0. 28l25E-4,0. 285E-7,0. 32625E-10  / 
DATA  (  CST(I) ,1=1,3  )  / 

. 5625E-4 , . 57E-7 , . 6512E-10  / 

IP  (  N.NE. 2. AND.N.GE. 1  )  GO  TO  3011 
IND=3 


l:-6 


3011 


301  i 

3031 

3041 

3051 

3061 


RETURN 
CONTINUE 
P=0 . 

IND-0 

A=0. 

KOM=0 . 

K=1 

IF  {  N.NE.l  )  GO  TO  3021 
W=X (2) -X (1) 

Z=Y(2)-Y{1) 

U  =X (3) -X (1) 

V  =Y (3) -Y (1) 

PSI1=V*W-U*Z 

IF  (  PSI1.GE.0.  )  GO  TO  3041 
P=-l. 

Tl=W 

W=U 

U=T3 

Tl=v 

r  z 

Z-=T1 

GO  TO  3041 
CONTINUE 
Y (N+l) =Y (1) 

X (N+l) =X (1) 

U  =X (2) -X (1) 

V  =Y (2) -Y (1) 

XK=X(1) 

YK=Y  (1) 

CONTINUE 
W=X (1) -X (N) 

Z=Y (1) -Y (N) 

CONTINUE 

DlSQ=W*W+Z*Z 

IF  (  D1SQ.GT.TAUSQ  )  GO  TO  3051 
IF  (  N.EQ.l  )  GO  TO  4011 
N=N-1 

IF  (  N . EQ . 2  )  RETURN 
GO  TO  3031 
CONTINUE 
D2SQ=U*U+V*V 

IF  (  D2SQ.GT.TAUSQ  )  GO  TO  3071 
IF  (  N.EQ.l  )  GO  TO  4011 
CONTINUE 
K=K+i 

U»X (F+l) -XK 
V*=Y  (K4-1)  -YK 
D2SQeU*U+V*V 
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IF  (  D2SQ.LE.TAUSQ  )  GO  TO  3061 
IF  (  K.EQ.(N-l)  )  RETURN 
3071  CONTINUE 

A=XK* (Y (K+l) -Y (N) ) 

BGD1~SQRT(2. *D1SQ) 

BGD2=SQRT(2.*D2SQ) 

3081  CONTINUE 

PSI1=V*W-U*Z 

CEE=U*W+V*Z 

AJO  =ATAN2(PSI1,CEE) 

KOM=KOM+AJO 

L=0. 

B*=.  5*  (X  (K)  *X  (K)  +Y  (K)  *Y  (K) ) 

IF  {  B.GT.APHl(IOP)  )  GO  TO  3111 
CAPG=0. 

3101  CONTINUE 

PI  =AJ0  /TWOPI-CAPG 
GO  TO  3621 
3111  CONTINUE 

G (1) * (W*X (K) +Z*Y (K) ) /BGD1 
G  (2)  =  (U*X (K) +V*Y  (K) ) /BGU2 
H(1)=(-Y(K)*W+X(K)*Z}/BGD1 
H (2) = (-Y (X) *U+X (X) *V) /BGD2 

IF  (  ABS(PSIl) .GT. <BGDl*BGn2*A3D8(IOP>) )  GO  TO  3241 
IF  (  CEE.LT.O.  )  GO  TO  3131 
IF  (  ABS(AJG) .LE.TAU  >  GO  TO  3121 
IF  {  G (1) .GE.O.  )  GO  TO  3121 
GO  TO  3241 
3121  CONTINUE 
Pl=0, 

GO  TO  3621 
3131  CONTINUE 

IF  (  ABS(PSIl) .LE.{.5*TAU*BGD1*BGD2)  )  IND»2 
IF  {  PSI1.LT.0.  )  GO  TO  3171 
P1  s.5*ERFCl(0fH(2)) 

GO  TO  3621 
3171  CONTINUE 

PI  . 5*ERFC1 (0#H (1) ) 

GO  TO  3621 
3^41  CONTINUE 

IF  (  B.LE. APH2 (IOP)  )  GO  TO  3301 
IF  {  G  (1)  .M.0.  )  GO  TO  3261 
IF  (  G{2) ,G2.0 .  )  GO  TO  3471 
G  (2)  *~G  (2) 

B(2)«-H(2) 

IF  {  ABvS{H<2))  .LS.A"H4(IGPj  )  GO  TO  3251 
L*.5*ERPC.1(G,-H(2) ) 

GO  TO  3461 


3251  CONTINUE 

L=.5+RTPII*H(2) 

GO  TO  3461 
3255  CONTINUE 

L=.5-RTPII*H(1) 

GO  TO  3461 
3261  CONTINUE 
G (1) =-G (1) 

H (1) =-H{l) 

IF  (  G ( 2 ) , LT . 0 .  )  GO  TO  3271 

IF  (  ABS(H(1) ) .LE.APH4 (IOP)  )  GO  TO  3255 

L= . 5*ERFC1 ( 0 ,H (1) ) 

GO  TO  3461 
3271  CONTINUE 
G (2)  =-G (2) 

H (2) =-H (2) 

IF  (  ABS (H(l) ) .LE.APH4 (IOP)  )  GO  TO  3291 
IF  (  ABS(H(2) ) .LE.APH4 (IOP)  )  GO  TO  3281 
L= . 5* (ERFC1 (0  #H (1) ) -ERFC1 (0  r H (2) ) ) 

GO  TO  3471 
3281  CONTINUE 

L=RTPII*H(2)-.5*ERF1{H(1) ) 

GO  TO  3471 
3291  CONTINUE 

IF  (  ABS ( H  ( 2 ) )  .  LE . APH4 ( IOP)  )  GO  TO  3295 
L« . 5*ERF1 (H (2) ) -RTPII*H (1) 

GO  TO  3471 
3295  CONTINUE 

L=RTPII* (H (2) -H(l) ) 

GO  TO  3471 
3301  CONTINUE 

CAPG=‘C1*<H'2)-H(1))-C2*{G(2)*H<2)-G(1)*H(1)) 
GO  TO  3101 
3461  CONTINUE 

IF  (  PSX1» LE.O*  )  GO  TO  3465 
L*L-1. 

AJOaPI+AJO 
GO  TO  3471 
3465  CONTINUE 

AJO“AJO  -PI 
3471  CONTINUE 

IF  (  B.GE.RSQ(IOP)  )  GO  TO  3501 

CAPE=AJ0 

CAPH».5*AJ0 

M«1 

F«0. 

A*?1“H{2)  -H  (1) 
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CIRCM=AJ1 

IF  (  I0P.EQ.3  )  GO  TO  3681 
IF  (  IOP.EQ.2  )  GO  TO  3701 
SUM=E (M) *AJ1 
3481  CONTINUE 
M=M+1 

H (2)  =H (2)  *G  (2) 

H  (1)  =H  (1)  *G  (1) 

T=H (2) -H (1) 

F=F+B 

CAPV= (F*CAPE+T) /M 
SUM=SUM+E (M) *CAPV 
IF {  M  .GE.  5  )  GO  TO  3491 

CAPE=CIRCM 
CIRCM=CAPV 
GO  TO  3481 
3491  CONTINUE 

PI  =L+EXP (- (B+ALNPI) ) * (CAPH-SUM) 
GO  TO  3621 
3501  CONTINUE 
Pl=L 

3621  CONTINUE 

IF  (  K.NE.N  )  GO  TO  3651 
IF  (  N.NE.l  )  GO  TO  3631 
P= ABS ( P+ABS { Pi ) ) 

RETURN 

3631  CONTINUE 
P=P-P1 

KOM=KOM/TWOPI 

A=.5*A 

IF  {  KOM.LT.O.  )  GO  TO  3641 
KO=INT {KOM+. 125  ) 

GO  TO  3645 
3641  CONTINUE 

KO*INT (KOM~ .125  ) 

3645  CONTINUE 

P-P4-FLOAT  (KO) 

RETURN 

3651  CONTINUE . 

W=U 

z«=v 

BGD1«BGD2 
XK*X (K+ 1) 

YK*Y(K+1) 

YKM1*Y (K) 

3661  CONTINUE 
K®K-fl 

U®X(K«fl)-XK 
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V=Y(K+1)-YK 

D2SQ=U*U+V*V 

IF  (  D2SQ.LE.TAUSQ  )  GO  TO  3661 
BGD2=SQRT(2.*D2SQ) 

P=P-P1 

A-A+XK* (Y (K+l) -YKM1  ) 

GO  TO  3081 
3681  CONTINUE 

SUM=E3 (M) *AJl 
3691  CONTINUE 
M=M+1 

H(2) -H{2) *G(2) 

H(1)«H(1)*G(1) 

T=H(2) -H(l) 

F=F+B 

CAPV= (F*CAPE+T) /M 
SUM=SUM+E3 (M) *CAPV 
IF  (  M.GE.15  )  GO  TO  3491 
CAPE-CIRCM 
CIRCM-CAPV 
GO  TO  3691 
3701  CONTINUE 

SUM=E2 (M) *AJ1 
3711  CONTINUE 
M=M+1 

H ( 2 ) =H ( 2 ) *G ( 2 ) 

H(1)«H{1)*G(1) 

T^H ( 2) -H (1) 

Ff=F+B 

CAPV-(F*CAPE+T)/M  , 

SUM«*SUM+E2  (M)  *CAPV 
IF  {  M.GE.10  )  GO  TO  3491 
CAPE-CIRCM 
CIRCM-CAPV 
GO  TO  3711 
4011  CONTINUE 
P-5. 

IND-1 

RETURN 

ENE 


FUNCTION  ERF1 (X) 

DIMENSION  A (4) ,B(4) ,P(8) ,Q(8) ,R(5) ,S(5) 

DATA  A/2. 426679 55230 53 2E0 2,  2 . 19792616182942E01 , 
1  6. 996 38348861914 E0 0,-3. 56 0984 37018 154E-2/ 

DATA  B/2 . 15058 8 7 586 986 1E0 2,  9 . 11649054045149E01 , 
1  1.50827976304078E01,  1. OOOOOOOOOOOOOOEOO/ 

DATA  P/3.00459261020162E02,  4. 51918953711873E02, 

1  3. 39320816734344 E0 2,  1. 52989285046940E02 , 

2  4. 3 16222 72220567 E01,  7 . 21175825088309E00 , 

3  5. 6419 5 517 4 789 74 E- 1,-1. 36 8 64 8 5738 2717E-7/ 
DATA  Q/3 . 00459260956983 E02,  7 . 90950925327898E02 , 

1  9.31354094850610E02,  6 . 38980264465631E02 , 

2  2. 7758 544474 3988E02,  7 . 70001529352295E01, 

3  1.278272731962 94E01,  1. OOOOOOOOOOOOOOEOO/ 
DATA  R/2 . 9 9610 70 7 70 3 54 2E-3 ,  4. 94730910623251E-2, 

1  2. 26 9 56 593539687 E-l,  2. 78661308609648E-1, 

0  9  9  VI  Q94ciQ7'}41  ft^F-9  / 

DATA  S/1 ! 0620 9230 52846 8E- 2,  1. 91308926107830E-1, 

1  1.05167510706793EOO,  1. 98733201817135E0O, 

2  1. OOOOOOOOOOOOOOEOO/ 

DATA  C/ 5. 64189583547756 E-l/ 


AX=ABS (X) 

X2=AX*AX 

IF  (AX.GE.0.5)  GO  TO  20 
TOP=A ( 4 ) 

BOT=B ( 4 ) 

DO  10  1=1,3 
J=4-I 

TOP=A{J) +X2*TOP 
10  BOT=B(J)+X2*BOT 
ERFl=X*TOP/BOT 
RETURN 
C 

20  IF  (AX.GT.4.0)  GO  TO  30 
TOP»P (8) 

BOT=Q (8) 

DO  21  1=1,7 
J*8-I 

TOP=P ( J) +AX*TOP 

21  BOT=Q { J) +AX*BOT 
ERFl=l.-EXP(-X2) *TOP/BOT 
IF  (X.LT.O.)  ERF1=-ERF1 
RETURN 

C 

30  ERFl=l . 

IF  (AX.GE. 5.54}  GO  TO  32 
TOP=R(l) 


BOT=S (1) 

DO  31  1=2,5 
TOP=R(I)+X2*TOP 

31  BOT=S (I) +X2*BOT 

ERF1C-TOP/  (X2*BOT) 
ERFl=l.-EXP (-X2) *ERF1/AX 

32  IF  (X.LT.O.)  ERFl=-ERFl 
RETURN 

END 


FUNCTION  ERFC1 { IND,X) 

DIMENSION  A(4) ,B(4) fP(8) ,Q(8) ,R(5) ,S(5) 

DATA  A/2. 426 6 79 5 5 230 53 2E0 2,  2.19792616182942E01, 
1  6 . 99638348861914E00 ,-3 . 56098437018154E-2/ 

DATA  B/2. 150 58 87586986 1E02,  9 . 11649054045149E01, 
1  1.50827976304078E01,  1. OOOOOOOOOOOOOOEOO/ 

DATA  P/3 . 004 59 26 102016 2E0 2 ,  4.51918953711873E02, 

1  3. 393208 16 734344 E0 2,  1. 52989285046940E02 , 

2  4. 31622272220 56 7 E01,  7 . 21175825088309E00 , 

3  5. 64195517478974E-1,-1. 36864857382717E-7/ 
DATA  Q/3.00459260956983E02,  7 . 90950925327898E02 , 

1  9.31354094850610E02,  6 . 38980264465631E02, 

2  2.77585444743988E02,  7 . 70001529352295E01, 

3  1.27827273196294E01,  1. OOOOOOOOOOOOOOEOO/ 

DATA  R/2. 99610707703542E-3 ,  4.94730910623251E-2, 

1  2.26956593539687E-1,  2.78661308609648E-1, 

2  2.23192459734185E-2/ 

DATA  S/1. 06 209 230 52846 8 E- 2,  1. 91308926107830E-1, 

1  1.O516751O706793EOO,  1.987332O1817135E00, 

2  1. OOOOOOOOOOOOOOEOO/ 

DATA  C/5.64189583547756E-1/ 


AX=ABS (X) 

X2=AX*AX 

IF  (AX. GE. 0.47)  GO  TO  20 
TOP=A(4) 

BOT=B{4) 

DO  10  1=1,3 
J=4-I 

TOP=A(J)4X2*TOP 
10  BOT=B(J)+X2*BOT 

ERFCl =1 . -X*TOP/BOT 
IF  (IND.NE.O)  ERFC1=EXP (X2) *ERFC1 
RETURN 
C 

20  IF  (AX.GT.4.0)  GO  TO  30 
TOP=P{8) 

BOT=Q(8) 

DO  21  1=1,7 
J=8-I 

TOP“P ( J) +AX*TOP 

21  BOT=Q(J)+AX*BOT 

ERFCl aTOP/BOT 
IF  (IND.EQ.0)  GO  TO  22 
IF  (X.LT.0.0)  ERFC1=2 . *EXP (X2) -ERFCl 
RETURN 

22  ERFClaEXP(-X2) *ERFC1 

IP  (X.LT.0.0)  ERFC1=2. -ERFCl 
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RETURN 


30  IP  (X.LE.-5.33)  GO  TO  32 
TOP=R(i) 

BOT=S (1) 

DO  31  1=2,5 
TOP=R(I)+X2*TOP 

31  BOT=S(I)+X2*BOT 

ERFC1= (C-TOP/ (X2*BOT) ) /AX 
IF  (IND.EQ.O)  GO  TO  22 
IP  (X.LT.0.0)  ERFC1=2 . *EXP (X2) -ERFCl 
RETURN 

32  ERFC1=2. 

IF  (IND.NE.O)  ERFC1=EXP(X2) *ERFC1 

RETURN 

END 


SUBROUTINE  SORT  III 
(FLOW  CHART  3,  page  43) 

Subroutine  SORT  III  Used  to  Eliminate  SDP  and/or  SCP  from  II 


CALL:  SORT  III  (x,  y,  N), 
where: 

x  is  the  array  of  abscissas  of  the  numbered  points  of  the  polygon  II.  The  array  is  dimensioned 
at  N.  Upon  return  to  the  calling  program,  P-2  (or  P-7),  the  array  of  abscissas  will  be 
reduced  by  the  number  of  consecutive  duplicate  points  SDP  and  SCP  eliminated.  The 
array  is  compacted.* 

y  is  the  array  of  ordinates  of  the  numbered  points  of  the  polygon  II.  The  array  is  dimensioned 
at  N.  Upon  return  to  the  calling  program,  P-2  (or  P-7),  the  array  of  ordinates  will  be 
reduced  by  the  number  of  points  deleted  due  to  SDP  and  SCP.  The  array  is  CU. 

N  is  the  number  of  points  initially  used  to  specify  the  polygon.  Upon  return  to  the  calling 
program,  P-2  (or  P-7),  N  will  be  reduced  by  the  number  of  points  that  were  eliminated. 


•Compact  hero  means  Uni  whenever  a  point  is  eliminated  all  subsequent  points  of  the  array  arc  moved  up  one 
location  in  the  array,  i.c.,  the  array  is  closed  up  (CU). 
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SUBROUTINE  SORT3  (  X,Y,N  ) 
DIMENSION  X(i),Y(l) 

DATA  CST/4. 9E-27  / 

3041  CONTINUE 

IF  (  N.LT.3  )  RETURN 

K=1 

L=2 

3051  CONTINUE 

U=X (1) -X (N) 

V=Y (1)  -Y (N) 

D2=U*U+V*V 

IF  (  D2.GT.CST  )  GO  TO  3061 
N=N-1 

IF  (  N.GT.2  )  GO  TO  3051 
RETURN 

3061  CONTINUE 

W=X  (L)  -X  (1) 

Z®Y (L) -Y (1) 

D1=W*W+Z*Z 

IF  (  Dl.GT.CST  )  GO  TO  3071 
L=L+1 

GO  TO  3061 
3071  CONTINUE 

IF  (  L. EQ . (K+l)  )  GO  TO  3091 

LM2=L-2 

N=N-LM2 

DO  3081  I»2,N 

Il<=LM2+I 

X(I)=X(I1) 

Y  (I)  =Y (II) 

3081  CONTINUE 
L=»2 

3091  CONTINUE 
T=V*W-U*Z 

SN» (4. *T*T) / (Dl*D2) 

IF  (  SN.GT.CST  )  GO  TO  3121 
3111  CONTINUE 
L=L+1 

IF  (  L.GT.N  )  GO  TO  3341 
3115  CONTINUE 

W“X (L) -X (1) 

Z=Y{L)-Y(1) 

D1=W*W+Z*Z 

IF  (  Dl.GT.CST  }  GO  TO  3091 
GO  TO  3111 
3121  CONTINUE 

IF  (  L . EQ . 2  )  GO  TO  3141 
LN2=L-2 
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N=N-LM2 

DO  3131  I=lrN 
Il=LM2+I 
X  (I)  =X  (II) 

Y  (I)  =Y  (II) 

3131  CONTINUE 

GO  TO  3041 
3141  CONTINUE 
K=2 
L=3 

GO  TO  3161 
3151  CONTINUE 
D1=D2 
W=U 

z=v 

3155  CONTINUE 
L=K+1 

3161  CONTINUE 

U=X(L)  -X(K) 

V=Y (L) -Y (K) 

D2=U*U+V*V 

IF  (  D2.GT.CST  )  GO  TO  3171 
3165  CONTINUE 
L=L+1 

IF  (  L.LE.N  )  GO  TO  3161 
N«K 

GO  TO  3251 
3171  CONTINUE 

IF  (  L.EQ. (K+l)  )  GO  TO  3191 
N=N-((L-1)-K) 

KPl=K+l 

I2=L-KP1 

DO  3181  X«KPl,N 
11=12+1 
X (I) ®X (II) 

Y(I)  BY (II) 

3181  CONTINUE 
L“KPl 

3191  CONTINUE 
T«V*W-U*Z 

SN“ (4 . *T*T) / (Dl*D2) 

IF  (  SN.GT.CST  )  GO  TO  3221 
3201  CONTINUE 
L°L+1 

IF  {  L.GT.N  )  GO  TO  3211 
U=X(L)-X(K) 

V=Y(L)-Y(K) 

D2=U*U+V*V 
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3211 

3221 

3231 

3241 

3251 

3255 

3261 


IF  (  D2.GT.CST  )  GO  TO  3191 
GO  TO  3201 
CONTINUE 
X  (K)  =X  (N) 

Y(K)=Y(N) 

N=K 

GO  TO  3251 
CONTINUE 

IF  (  L.EQ. (K+l)  )  GO  TO  3241 

I2=L-1-K 

N=N-I2 

LM2=L-2 

DO  3231  I=K,N 
11=12+1 
X(I)=X(I1) 
y(I)=Y(Il) 

CONTINUE 
W=X (K) -X (K-l) 

Z=Y(K)-Y(K-1) 

D1=W*W+Z*Z 

IF  (  Dl.GT.CST  )  GO  TO  3155 
K=K-1 

IF  (  K.LT.2  )  GO  TO  3041 
W=X  (K)  -X  (K-l) 

Z=Y (K) ~Y (K-l) 

Dl»W*W+Z*Z 

L-K+l 

GO  TO  3165 

CONTINUE 

K»K+1 

IF  (  K.LT.N  )  GO  TO  3151 

GO  TO  3255 

CONTINUE 

U»X(N)-X(N-1) 

V=Y(N)-Y(N-1) 

D2»U*U+V*V 

IF  {  D2.LE.CST  )  GO  TO  3261 

CONTINUE 

W*XU)-X(N) 

Z*Y (1)  -Y (N) 

D1-W*W+Z*Z 

IF  (  Dl.Lfi.CST  )  GO  TO  3261 
T=V*W-U*2 

SN« (4 . *T*T) / (D1*D2) 

IF  {  SN.GT.CST  )  GO  TO  3351 

CONTINUE 

N=N-1 

IF  (  N.GT.2  )  GO  TO  3251 
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RETURN 

3341  CONTINUE 
N=2 

RETURN 

3351  CONTINUE 
D2=D1 
U=W 
v=z 

W=X (2) -X (1) 

Z=Y (2) -Y (1) 

Dl=W*W+Z*Z 

T=V*W-U*Z 

SN=(4,*T*T)/(D1*D2) 

IF  (  SN.GT.CST  )  RETURN 
L=3 

GO  TO  3115 
END 


SUBROUTINE  VALR-7 
(FLOW  CHART  4,  page  44) 

Subroutine  VALR-7  Used  to  Compute  p(S),  where  S  has  no  SAR(s),  (See  page  31) 


CALL:  VALR-7  (x,  y,  M,  p,  IOP,  a,  IND),* 
where: 

x  is  the  input  array  of  abscissas  for  S.  Dimensioned  at  M  +  1. 

y  is  the  input  array  of  ordinates  for  S.  Dimensioned  at  M  +  1. 

M  is  the  number  of  input  points  for  S.  When  M  =  1,  IBND  over  an  angular  region  is  com¬ 

puted.  Three  input  points  in  counterclockwise  order  are  used  to  specify  the  region  with 
the  vertex  at  (1). 

p  is  the  location  where  the  function  value  for  p(S)  will  be  returned.''' 

IOP  is  set  by  the  user  to  1,  2,  or  3  for  approximately  3,  6,  or  9-decjmal-digit  accuracy, 

respectively,  in  p(S). 

a  is  the  location  where  the  value  of  the  function  a(S)  is  returned.  The  absolute  value  of 
a  gives  the  area  of  S. 

IND  is  an  error  indicator  normally  set  to  zero.  If  PAR(s)  are  detected  by  VALR-7,  then 
IND  is  set  to  two  and  the  result  for  p(S)  is  most  likely  wrong,  unless  M  =  1.  See  Flow 
Chart  4-24,  20,  21,  22.  VALR-7  should  never  be  used  alone  if  SAR(s)  are  a  possibility, 
unless  M  =  1.  IfM  =  2orM<l,  then  IND  «  3  and  an  EXIT  is  made.  Such  M  are  not 
allowed. 

This  routine  requires  computation  of  erf  (x)  and  erfc  (x)  which  are  defined  on  pages  5,  28  and 
29.  We  have 


ERF  1  (x)  =  erf(x),  ERFC  1  (0,  x)  =  erfc(x), 

where  the  subroutine  listings  for  these  functions  are  given  on  pages  F-12  to  F-15.  They  are  identical 
to  the  NSWC(DL)  math  library  functions  ERF  and  ERFC  as  of  June  1980  which  are  based  on  the 
reference  below. 

Cody,  W.  J.,  Rational  Chebyshev  Approximations  for  the  Error  Function,  Mathematics  of 
Computation,  v.  23  (1969),  pp.  631-637. 


*  vVe  use  p,  a,  M  here  In  place  of  P,  A,  N  to  avoid  ambiguity  with  results  in  P-7,  if  SORT  1  is  used  with  VALR-7. 
tThe  IBND  over  S,  p(S),  will  be  positive  if  S  is  PO  and  it  will  be  negative  if  S  is  NO. 
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SUBROUTINE  VALR7  (  X ,Y ,N, P, IOP, A, IND  ) 
DIMENSION  RSQ ( 4 ) 

DIMENSION  X (1) ,Y<1) ,G(2) ,H(2) 

DIMENSION  £(5) ,E2{10) ,E3{15) 

DIMENSION  APH1 (3) ,APH2(3) ,APH4(3) ,CST(3) 
REAL  L 


DATA  TWOPI/6 .2831  85307  17958  / 
DATA  ALNPI/1. 1447  29885  84940  / 
DATA  Cl/. 28209  47917  73877  / 
DATA  C2/. 15915  49430  91895  / 
DATA  TAU/7.E-14  / 

DATA  (  E(I),I=1,  5)  / 

. 88 5777518 5728 95E+00  , 

. 759305502082485E+00  , 

. 69523 20 92435207 E-01  / 
DATA  (E2 (I)  , 1=1,  10)  / 

.8 862264 700166 3 2E+00  , 

. 885348820003892E+00  , 

. 421821 197 160099E+00  , 

. 905057384150449E-01  , 

. 430895168984138E-02  , 
DATA  (E3(I),I=1,  15)  / 

.886226924931465E+00  , 

. 8862237 33186722E+00  , 
.442851899328569E+00  , 
.145060043403014E+00  , 

. 309199295521210E-01  , 

. 32 494 45431 71185E-02  , 

. 105787574480633E-03  f 
. 408335517232165E-06  / 
DATA  (  APHl(I) f 1=1,3  )  / 
2.02E-7,2.08E-13,2.71E-19  / 
DATA  (  APH2(I) ,1=1,3  )  / 
1.22E~2,1.23E-4,l,34E-6  / 


- . 981151952778050E+00 
353644980686977E+00 


- . 999950714561036E+00 
- . 66 06112390 433 57E+00 
- . 222898055667208E+00 
-.254906111884287E-01 
-*  .32337  7239G93247E-03 

- .999999899776252E+00 
- . 666626670510907E+00 
- . 265638206366025E+00 
-» 71490 983 7799889E-01 
- . 112323532148441E-01 
- . 704260243309096E-03 
-.971864864160461E-05 


DATA  (  APH4(I) ,1=1,3  )  / 

.6962E-1,  .6990E-2,  .7311E-3  / 

DATA  RTPII/. 56418  95835  4776  / 

DATA  (  RSQ  (I), 1=1, 3  )  / 
6.0516,12.60605  ,19.201924  / 

DATA  (  CST(I) , 1=1, 3  )  / 

. 5625E-4 , . 57E-7 , . 6512E-10  / 

IP  (  N.NE.2. AND.N.GE. 1  )  GO  TO  3061 

IND=3 

RETURN 

CONTINUE 


P=0. 


IND=0 


IP  (  N.NE.l  )  GO  TO  3071 


/ 
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K=1 

A=0. 

W=X (2) -X (1) 

Z=Y (2) -Y (1) 

U  =X (3) -X (1) 

V  =Y(3)-Y(1) 

PSI1=V*W-U*Z 

IP  (  PSI1.GE.0.  )  GO  TO  3081 
P=-l. 

Tl=W 

W=U 

U=Tl 

Tl-V 

V=Z 

Z=Tl 

GO  TO  3081 
3071  CONTINUE 

CALL  SMP7  (  N,A,X,Y  ) 

IF  (  ABS (A  } .LE.CST(IOP)  )  RETURN 
K=1 

W=X(1)-X(N) 

Z=Y(1)-Y(N) 

U  =X (2) -X (1) 

V  =Y (2) -Y (1) 

X (N+l) =X (1) 

Y (N+l) =Y (1) 

3081  CONTINUE 

BGD1 =SQRT {  2, * (W*W+Z*Z) ) 

0GD2 =SQRT (  2.* (U*U+V*V) ) 

3091  CONTINUE 
L»0, 

B=.5*(X(K)*X(K}+Y(K)*Y(K)) 

IP  (  B.GT. APHl ( IOP)  )  GO  TO  3111 
CAPG-0, 

3101  CONTINUE 

T1*V*W-U*Z 
T2"U*W+V*Z 
PHIK*ATAN2 (Tl ,T2) 

PI  -PHIK/TWOPI-CAPG 
GO  TO  3621 
3111  CONTINUE 

G  (1)  “  (W#X  (K)  +3*Y  (K) )  /BGDl 
G (2) - (U*X <K) +V*Y (K) )/BGD2 
H (1) ® (-Y (K) *W+X (K) *Z) /BGDl 
H(2)«(-Y(K)*U+X(K)*V)  /BGD2 
SN- (2. * (V*W-U*Z) } / (BGDl*BGD2) 

IF  (  ABS(SN) .GT.CST(IOP)  )  GO  TO  3241 
CN*G{1)  *G(2)+H{1)  *H(2) 
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IF  (  CN.LT.O.  )  GO  TO  3131 
IF  (  ABS(SN) .LE.TAU  )  GO  TO  3121 
IF  (  G (1) .GE. 0.  )  GO  TO  3121 
GO  TO  3241 
3121  CONTINUE 
P1=0. 

GO  TO  3621 
3131  CONTINUE 

IF  (  ABS(SN) .LE.TAU  )  IND=2 
IF  (  SN  . LT . 0 .  )  GO  TO  3171 
PI  =.5*ERFC1{0,H(2) ) 

GO  TO  3621 
3171  CONTINUE 

Pi  =-. 5*ERFC1 (0 ,H (1) ) 

GO  TO  3621 
3241  CONTINUE 

IF  (  B.LE. APH2 (IOP)  )  GO  TO  3301 
SN=B*SN 

IF  (  G ( 1 ) . LT . 0 .  )  GO  TO  3261 
IF  (  G (2) .GE.O.  )  GO  TO  3471 
G (2)  =~G (2) 

H  (2)  s=-H  (2) 

IF  (  ABS ( H ( 2 ) ) . LE . APH4 ( IOP)  )  GO  TO  3251 
L=.5*ERFCl (0,-H{2) ) 

GO  TO  3461 
3251  CONTINUE 

L= . 5+RTPII*H (2) 

GO  TO  3461 
3255  CONTINUE 

L= . 5-RTPI I *H ( 1 ) 

GO  TO  3461 
3261  CONTINUE 
G  (1)  «-G  (1) 

H(l)=-H(l) 

IF  (  G (2) . LT . 0 .  )  GO  TO  3271 

IF  (  ABS  (H (1) )  . LE.  APH4  (IOP)  )  GO  TO  3255 

L=.5*ERFC1(0,H(1)) 

GO  TO  3461 
3271  CONTINUE 
G  (2)  b~G  (2) 

H  (2)  =-H (2) 

IF  (  ABS ( H { 1 } ) . LE . APH4 { IOP)  )  GO  TO  3291 
IF  (  ABS (H (2) } . LE* APH4 ( IOP)  )  GO  TO  3281 
L“. 5* (ERFCl (0,H(1) ) -ERFC1 (0, H(2) ) ) 

GO  TO  3471 
3281  CONTINUE 

L=RTPII*H(2)-.5*ERF1(H(1)) 

GO  TO  3471 
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3291  CONTINUE 

IF  (  ABS(H(2)) .LE.APH4 (IOP)  )  GO  TO  3295 
L=. 5*ERF1 (H (2) ) -RTPII*H (1) 

GO  TO  3471 
3295  CONTINUE 

L=RTPII*(H(2)-H(1) ) 

GO  TO  3471 
3301  CONTINUE 

CAPG=C1* (H ( 2 ) — H ( 1 ) ) -C2* (G(2)*H(2)-G(1)*H(1)) 
GO  TO  3101 
3461  CONTINUE 
SN=-SN 

IF  (  SN.LE.O.  )  GO  TO  3471 
L=L-1. 

3471  CONTINUE 

IF  (  B.GE.RSQ(IOP)  )  GO  TO  3501 
CN=G (1) *G (2) +H (1) *H(2) 

AJO  =ATAN  2 ( SN , CN ) 

CAPE=AJ0 
CAPH=. 5* AJO 
M=1 
F=0. 

AJl-H (2) -H (1) 

CIRCM=AJ1 

IF  (  IOP.EQ.3  )  GO  TO  3681 
IF  (  IOP. EQ. 2  )  GO  TO  3701 
SUM=E(M) *AJ1 
3481  CONTINUE 
M=M+1 

H (2)  =H(2) *G (2) 

H  (1)  =H(1)  *G(1) 

T=H{2)-H(1) 

F=F+B 

CAPV={F*CAPE+T)/M 
SUM*SUH+E(M) *CAPV 
IF (  M  .GE.  5  )  GO  TO  3491 

CAPE*CIRCM 
CIRCM»CAPV 
GO  TO  3481 
3491  CONTINUE 

PI  “L+EXP{- (B+ALNPI) ) * (CAPH-SUM) 

GO  TO  3621 
3501  CONTINUE 
P1*L 

3621  CONTINUE 

IF  (  K.NE.N  )  GO  TO  3651 
IF  (  N.NE.l  )  GO  TO  3631 
P® ABS ( P+ABS ( Pi ) ) 


RETURN 

3631  CONTINUE 
P=P-P1 

IF  (  A.LT.O.  }  GO  TO  3641 
P=P+1. 

RETURN 

3641  CONTINUE 
P=P-1. 

RETURN 

3651  CONTINUE 
K=K+1 
W=U 
Z=V 

U=X (K+l) -X (K) 

V=Y (K+l) -Y (K) 

BGD1=BGD2 

BGD2=SQRT (  2 . * (U*U+V*V) ) 
P=P-P1 
GO  TO  3091 
3681  CONTINUE 

SUM=E3 (M ) *AJ1 
3691  CONTINUE 
M=M+1 

H ( 2 ) =H ( 2 ) *G ( 2 ) 

H (1) =H(1) *G (1) 

T=H (2) -H (1) 

F=F+B 

CAPV=(F*CAPE+T)/M 
SUM=SUM+E3 (M) *CAPV 
IF  (  M.GE.15  )  GO  TO  3491 
CAPE=CIRCM 
CIRCM=CAPV 
GO  TO  3691 
3701  CONTINUE 

SUM=E2 (M) *AJ1 
3711  CONTINUE 
M*>M+1 

H(2)»H(2)*G(2) 

H(l)  =H(1)  *G(1) 

T“H(2) -H  (1) 

F*F+B 

CAPV® (F*CAPE+T) /M 

SUM=SUM+E2  <M) *CAPV 

IF  (  M.GE. 10  )  GO  TO  3491 

CAPE-CIRCM 

CIRCH®CAPV 

GO  TO  3711 


MASTER  SUBROUTINE  P-7 
(FLOW  CHART  5,  page  A-17) 

SUBROUTINE  P-7  is  Used  for  Computing  P(TI)  *or  an  Arbitrary  Polygon  II* 


CALL:  P-7  (x,  y,  N,  P,  ICV,  IND,  IOP,  A), 
where: 

x  is  the  array  of  abscissas  of  the  numbered  points  of  II.  x  is  dimensioned  at  N  +  1 . 

y  is  the  array  of  ordinates  of  the  numbered  points  of  >1.  y  is  dimensioned  at  N  +  1 . 

N  is  the  number  of  points  specifying  II,  except  if  N  =  1  when  the  IBND  over  an  angular 

region  is  computed.  Three  input  poin-.  are  needed,  for  N  =  1,  given  in  counterclock¬ 
wise  order,  with  the  vertex  at  point  one. 

P  is  the  location  where  the  value  of  P(II)  is  returned. 

ICV  must  be  specified  as  an  integer  by  the  user  according  to  the  list  below: 

ICV  =  0  II  is  simple  or  of  S  type  with  no  SAR(s)  (pages  1 2, 3 1).  VALR-7  used  alone. 

ICV  =  1  II  is  in  (S).  SORT  W  used  with  VALR-7. 

ICV  =  2  II  is  in  {II}.  SORT  I  is  used  to  search  for  duplicate  points**  of  n  in  in¬ 

creasing  digital  oruer  from  point  (2)  to  point  (N).  II  is  numbered  with  the 
Qf-option  (see  page  14),  so  II  is  decomposed  into  simple  polygons,  Sl, 

S2 . SJ.  SORT  II  is  not  needed.  VALR-7  is  used  to  find  p(S'),  which  are 

summed  in  SORT  I  to  g  vo  P(I1). 

ICV  =  -2  IT  is  in  (11).  SORT  I  is  used  to  search  for  duplicate  points  of  II  in  decreasing 
digital  order  from  point  (N-!)  to  point  (1).  11  is  numbered  with  the 
a-option. 

ICV >3  I?  is  in  'll).  SORT  I  is  used  to  search  for  duplicate  points  in  increasing 
digit  order  of  the  numoored  points  from  point  (2)  to  point  (N).  II  is 
n  nbered  with  the  0-option  (see  page  24),  so  II  is  decomposed  into  S  type 
elements  S1 ,  ...,  SL.  These  elements  require  SORT  II  to  eliminate  any  SCP. 
so  that  VALR  7  can  be  used  on  cadi  S*  to  obtain  p(S*),  whidt  are  summed 
in  SORT  i  to  give  Pfil). 

ICV  <  0  Tins  has  the  same  function  as  ICV  =  3,  except  that  SORT  1  scardies  for 
^-2  duplicate  points  of  II  in  decreasing  digital  order  of  the  numbered  points 
starting  at  point  (N  -  I )  and  finishing  at  point  ( I ). 


‘See  footnote  I ,  page  I  for  definition  of  an  arbitrary  polygon. 

“Duplicate  points  arc  not  to  be  confused  with  SDPIs),  see  pages  43  and  A-S. 


F-27 


Generally  ICV  =  3  is  preferable  to  ICV  =  2  and  ICV  =  -3  is  preferable  to  ICV  =  -2, 
because  the  computing  time  may  be  less  since  often  fewer  angular  regions  of  II  will 
need  processing. 

IND  is  an  error  indicator.  It  is  normally  set  at  zero.  However,  if  VALR-7  is  used  alone 
(ICV  =  0)  on  a  polygon  containing  PAR(s),  then  IND  is  set  to  two  and,  unless  N  =  1 , 
the  result  for  P  is  probably  wrong.  This  will  never  occur  if  SORT  III  or  SORT  I  and 
SORT  II  are  used  to  eliminate  SAR(s)  before  using  VALR-7,  provided  II  is  in  {S}. 
If  N  is  not  set  to  one  or  greater  than  two,  as  an  integer,  then  IND  is  set  to  three  with 
direct  exit  from  VALR-7.  Such  N  are  not  allowed. 

IOP  is  an  accuracy  parameter.  It  is  set  by  the  user  to  1,  2,  or  3  for  approximately  3,  6,  or 
9-decimal  digits  of  accuracy  in  P(II). 

A  is  the  location  where  A(II)  is  returned.  |A|  gives  the  area  of  II.  (See  Appendix  D, 
also  (46).) 


SUBROUTINE  P7  (  X,Y,NB,P,ICV,IND,IOP,A  ) 
DIMENSION  X{1)  ,Y(1) 

IF  (  NB.NE. 2. AND.NB.GE.l  )  GO  TO  3031 

IND=3 

RETURN 

3031  CONTINUE 
N=NB 

IF  (  N  . EQ.l  )  GO  TO  3041 

IF  (  ICV  . EQ.O  )  GO  TO  3041 

IF  (  ICV, EQ.l  )  GO  TO  3061 

CALL  SORT!  (  X,Y,N,P, ICV, IND, IOP,A  ) 

RETURN 

3041  CONTINUE 

CALL  VALR7  (  X,Y,N,P, IOP, A, IND  ) 

RETURN 

3061  CONTINUE 

CALL  SORT 3  (  X,Y,N  ) 

IF  (  N  .GT.2  )  GO  TO  3071 
A-Q. 

IND=0 

P=0. 

RETURN 

3071  CONTINUE 

CALL  VALR7  (  X,Y,N,P,IOP,A,IND  ) 

RETURN 

END 


SUBROUTINE  SORT  I  (See  Appendix  A) 

(FLOW  CHART  6,  page  A-18) 

Subroutine  SORT  I  Used  to  Decompose  II  Into  S  or  S  Type  Elements 

CALL:  SORT  I  (x,  y,  N,  P,  ICV,  IND,  IOP,  A), 
where : 

x  is  the  array  of  abscissas  of  the  numbered  points  specifying  the  polygon,  IT.  x  is  dimen¬ 
sioned  at  N  +  1 . 

y  is  the  array  of  ordinates  of  the  numbered  points  specifying  the  polygon,  II.  y  is  dimen¬ 
sioned  at  N  +  1 . 

N  is  the  number  of  points  numbered  on  the  polygon. 

P  is  the  location  where  P(IT)  is  returned. 

ICV  is  a  user  specified  integer  according  to  the  listing  below: 

ICV  =  2  for  a  polygonal  element  of  the  class  {11}  (see  page  1).  SORT  I  searches  for 
duplicate  points  of  II  in  increasing  digital  order  from  point  (2)  to  point  (N). 
n  is  specified  by  numbering  points  of  II  according  to  the  a-option.  (See 
page  14.)  II  is  decomposed  into  S1,  ....  SJ.  VALR-7  is  called  to  compute 
p(S').  These  quantities  are  summed  to  give  P(II). 

ICV  =  -2  for  a  polygonal  element  of  (11).  SORT  I  searches  for  duplicate  points  of  II 
in  decreasing  digital  order  from  point  (N  -  1)  to  point  ( I).  II  is  numbered 
according  to  the  a-option. 

ICV  >3  for  a  polygonal  element  of  {II}.  SORT  I  proceeds  in  the  same  way  as  for 
ICV  =  2,  except  that  II  is  numbered  according  to  the  0-option  rather  than 
the  a-option  (see  page  24).  The  0-option  numbering  requires  that  SOR  T  II 
be  used  to  eliminate  SCP  in  any  of  the  S  elements  obtained  from  the 
decomposition  of  II  by  SORT  I  (sec  Flow  Chart  6.  page  A-18).  11  is  de¬ 
composed  after  using  SORT  II  into  Sl,  ....  ST  VALR-7  is  called  to  com¬ 
pute  each  p(Si).  The  p(S')  are  summed  to  give  P(fl). 

ICV  <0  for  a  polygonal  element  of  (II).  SORT  I  proceeds  in  the  same  way  as  for 

^  -2  ICV  =  -2  except  that  II  is  numbered  according  to  the  0-option  rather  than 

the  a-option.  The  0-option  numbering  requires  that  SORT  11  be  used  to 
eliminate  SCP  in  any  of  the  S  elements  obtained  from  the  decomposition 
of  II  by  SORT  I. 

If  N  docs  not  differ  using  the  a  or  0-option,  then  ICV  *  ±2  is  preferable  to  ICV  =£  t2.  How¬ 
ever,  if  N  is  reduced  by  using  the  0-option,  then  ICV  =*  ±2  is  preferable  since  fewer  calls  to  VALR-7 
will  be  needed. 
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SUBROUTINE  SORT  I  (Continued) 


IND  is  an  error  indicator.  Normally  it  is  set  to  zero.  If  a  rr-angular  region,  PAR,  is  detected 
by  VALR-7,  IND  is  set  to  two,  and  p(S)  is  very  likely  wrong,  unless  N  =  1 ;  consequently 
also  P(II)  will  be  wrong.  The  routine  P-7  is  designed,  if  properly  used,  so  that  this 
cannot  happen  under  the  a-option,  nor  can  it  occur  under  the  ^-option  since  SORT  II 
removes  SCP  before  a  call  is  made  to  VALR-7  (see  Flow  Chart  6).  If  N  =£  1  or  is  not 
greater  than  two,  as  an  integer,  IND  is  set  to  three  and  an  exit  is  made.  Such  values 
of  N  are  not.  allowed. 

IOP  is  set  by  the  user  to  l ,  2,  or  3  to  obtain  approximately  3,  6,  or  9-decimal-digit  accuracy, 
respectively,  for  P(fl). 

A  is  the  location  where  the  A- function  value  for  n  is  returned.  The  area  of  II  is  given 
by  i A !  (see  SMP-7,  pages  9,  F-37). 


SUBROUTINE  SORTl  (  X,Y,N,P,ICV 
DIMENSION  X(1),Y(1) 

DATA  CST/5.E-14  / 

P=0. 

A=0. 

IC-IABS ( ICV) 

IF  (  ABS (X  (N  )-X  (1) ) .GT.CST  ) 
IF  (  ABS (Y  (N  ) -Y  (1)). GT.CST  ) 
GO  TO  2321 
2311  CONTINUE 
N=N+1 

2321  CONTINUE 
X{N)=X(1) 

Y(N)*Y(1) 

JlST=2 

11=2 

2331  CONTINUE 

IF  (  ICV.GT.O  )  GO  TO  2361 
NUMP1=N+1 

DO  2351  J1=J1ST,N 
J  =NUMPl-Jl 
JP1=J+1 

DO  2341  K=JPl,N 

IF  (  ABS (X  (J  )-X  (K  )). GT.CST  ) 

IF  (  ABS (Y  (J  ) -Y  (K  )). GT.CST 

IST=J 

IEN=K 

J1ST=N-K+1 

IF  (  K.EQ.N  )  J1ST=2 
LST=IST+1 
GO  TO  2531 
2341  CONTINUE 
2351  CONTINUE 
2361  CONTINUE 

DO  2521  1=11, N 

IM1=I-1 

DO  2511  Kl=*l, IMl 
K=I-K1 

IF  (  ABS (X  (I  )~X  (K  )). GT.CST 

IF  (  ABS (Y  {I  } -Y  (K  )). GT.CST 

IST=K 

IEN=I 

I1=K 

LST=IST 

IF  (  K.NE.l  )  GO  TO  2531 
11  =  2 

LST=LST+1 
GO  TO  2531 
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IND, IOP, A  ) 


GO  TO  2311 
GO  TO  2311 


GO  TO  2341 
GO  TO  2341 


GO  TO  2511 
GO  TO  2511 


2511  CONTINUE 
2521  CONTINUE 
2531  CONTINUE 

NUM1=IEN-IST 

NSAV=NUM1 

IF  (  NUM1.LE.2  )  GO  TO  2575 

IF  (  IC.EQ.2  )  GO  TO  2563 

CALL  SORT2  {  X  (1ST) ,Y  (1ST) ,NUM1  ) 

IF  (  NUM1.LT.3  )  GO  TO  2575 
2565  CONTINUE 

CALL  VALR7  (  X  (1ST)  #Y  (1ST)  ,NUM1#SMP,  IOP,SMA,IND  ) 

IF  (  IND.EQ.2  )  RETURN 

A=A+SMA 

P=P+SMP 

X  (1ST)  =X  (IEN) 

Y ( 1ST) =Y (IEN) 

2575  CONTINUE 

IF  (  IEN.NE.N  )  GO  TO  2577 
IF  (  IST.EQ.l  )  RETURN 
X(IST)=X(N) 

Y  (1ST)  SY  (N) 

N=IST 
GO  TO  2331 
2577  CONTINUE 
N=N-NSAV 

DO  2581  L=LST,N 
K*L+NSAV 
X  (L)  °X  (K) 

Y (L)  «Y(K) 

2581  CONTINUE 

GO  TO  2331 
END 
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SUBROUTINE  SORT  fl  (See  Appendix  A) 

(FLOW  CHART  7,  page  A- 19) 

Subroutine  SORT  II  Used  to  Eliminate  Successive  Colinear  Points  in  S 

CALL:  SORT  II  (x,  y,  M), 
where: 

x  is  the  array  of  abscissas  of  the  numbered  points  of  the  polygon  S.  The  array  is  dimensioned 
at  M.  Upon  return  to  the  calling  program  SORT  1,  the  array  of  abscissas  will  be  reduced 
by  the  number  of  points  deleted,  because  of  SCP.  The  x  array  is  compacted  or  closed  up. 

y  is  the  array  of  ordinates  of  the  numbered  points  of  the  polygon  S.  The  array  is  dimen¬ 
sioned  at  M.  Upon  return  to  the  calling  program  SORT  I,  the  array  of  ordinates  will  be 
reduced  by  the  number  of  points  deleted,  because  of  SCP.  The  reduced  y  array  is 
compacted  or  dosed  up. 

M  is  the  number  of  points  of  the  polygon  S  that  are  numbered.  Upon  return  to  the  calling 
program.  SORT  I,  M  will  be  reduced  by  the  number  of  successive  colinear  points  that 
were  eliminated 


SUBROUTINE  SORT2  (  X,Y,N  ) 
DIMENSION  X(1),Y(1) 

DATA  CST/4.9E-27  / 

K=1 

L=2 

U=X (1) -X (N) 

V-Y (1) -Y (N) 

D2=U*U+V*V 
3051  CONTINUE 

W=X (L) -X (1) 

Z=Y (L) -Y (1) 

Dl=W*W+Z*Z 

T=V*W-U*Z 

SN=(4.*T*T)/(D1*D2) 

IF  (  SN.GT.CST  )  GO  TO  3071 
L=L+1 

IF  (  L.LT.N  )  GO  TO  3051 
N=2 

RETURN 

3071  CONTINUE 
K=2 

iF  (  L.NE.2  )  GO  TO  3081 
L=3 

GO  TO  31U 
3081  CONTINUE 
LM2=L-2 
N=N- (LM2) 

DO  3091  1*1, N 

I1=LM2+I 

X(I)»=X(I1) 

Y ( I) =Y ( II) 

3091  CONTINUE 
3101  CONTINUE 
L«K+1 

W^X (K) -X (K~l) 

Z=Y (K) -Y (K-l) 

Dl=W*W+Z*Z 
3111  CONTINUE 

U=X(L)-X{K) 

V=Y (L) -Y (K) 

D2=U*U+V*V 

T=V*W-U*Z 

SN° (4. *T*T) / (D1*D2) 

IF  (  SN.GT.CST  )  GO  TO  3121 
L=L+1 

IF  (  L.LE.N  )  GO  TO  3111 
X  (K)  =X  (N) 

Y{K)=Y(N) 
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N=K 

GO  TO  3151 
3121  CONTINUE 

IF  (  L.EQ. (K+l)  )  GO  TO  3171 

LM2-L-2 

I3=LM2- (K-l) 

N=N-I3 

DO  3131  I=Kf N 
11=13+1 
X(I)=X(I1) 
y (i) =y (ii) 

3131  CONTINUE 
K=K+1 

IF  (  K.LT.N  )  GO  TO  3101 
3151  CONTINUE 

U=X (N) -X (N-l) 

V=Y (N) -Y (N-1) 

D2=U*U+V*V 
3161  CONTINUE 

W=X (1) -X (N) 

Z=Y(1)-Y(N) 

Dl=W*W+Z*Z 

T=V*W-U*Z 

SN=(4.*T*T)/{D1*D2) 

IF  (  3N.LE.CST  }  GO  TO  3165 
RETURN 

3165  CONTINUE 
N=N~1 
RETURN 

3171  CONTINUE 
K=K+1 

IF  (  K.GE.N  )  GO  TO  3161 

D1*D2 

W'=U 

z=v 

L=K+1 

GO  TO  3111 
END 


SUBROUTINE  SMP-7 
(No  flow  chart  given) 

SMP-7  is  Used  to  Compute  the  a-Function* 

CALL:  SMP-7  (M,  a,  x,  y)t, 
where: 

M  is  the  number  of  input  points  specifying  the  polygon, 
a  is  the  location  to  which  the  a-function  is  returned, 
x  is  the  array  of  input  abscissas.  Dimensioned  at  M. 
y  is  the  array  of  input  ordinates.  Dimensioned  at  M. 

(See  Appendix  D  for  value  of  a  in  the  wz-plane.) 


‘The  expression  used  to  compute  the  a-function  is  given  by 
i  N 

a  *i(y^i  y<)  °  Xn-  yN*t  “  y|*  (Sec  Appendix  D) 

"  I 

(The  area  of  the  input  polygon  is  given  by  la  j.) 

I'M,  a  are  used  in  place  of  N,  A  to  ovoid  confusion  with  the  latter  quantities  in  P-7  when  it  calls  SORT  t,  and 
SORT  I  in  turn  calls  VALR-7.  See  How  Giarts  S  and  6,  pages  (A- 1 7. 18). 
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SUBROUTINE  SMP7  (  NB,ANS,X,Y  ) 

DIMENSION  X{1)  ,Y(1) 

IF  (  NB.GT.3  )  GO  TO  3151 

ANS=. 5*  ((X{2)-X{1)}*(Y(3)-Y(1))-{X(3)-X(1))*(Y(2)-Y(1))) 
RETURN 

3151  CONTINUE 
NBM1=NB-1 

ANS=X (1) * (Y (2) -Y (NB)  ) +X (NB) * (Y (1) -Y (NBM1) ) 

DO  3161  1=2, NBM1 

ANS=ANS+X(I)  *(Y(I+1)-Y(M)  ) 

3161  CONTINUE 

ANS=. 5*ANS 

RETURN 

END 


APPENDIX  G 

TRIANGLE  CHECKOUT  PROGRAM  WITH  DREZNER 
(No  Flow  Chart) 


SUBROUTINE  DZ 

TRIANGLE  CHECKOUT  PROGRAM  with  DREZNER  (See  page  47) 

CALL:  DZ  (x,  y,  N,  P,  A), 
where : 

x  is  the  array  of  abscissas  of  the  points  specifying  polygon  II.  x  is  dimensioned  at  N. 
y  is  the  array  of  ordinates  of  the  points  specifying  II.  y  is  dimensioned  at  N. 

N  is  the  number  of  points  specifying  II. 

P  is  the  location  where  P(FI)  is  returned. 

A  is  the  location  where  A(I1)  is  returned. 

This  subroutine  decomposes  II  into  N-2  triangles  AJ  with  the  vertices  given  by  (1),  (j), 
(j  +  1),  j  =  2,  ....  N  -  1 .  P(Aj)  is  computed  by  DZ  -  1  and  A(Aj)  by  SMP-7;  the  results  are  summed 
in  DZ,  i.e.,  P(n>  =  P(Aj),  A(II)  =  Sf  A(Aj). 

This  routine  requires  computation  of  erf  (x)  and  erfc  (x)  which  are  defined  on  pages  5,  28  and 
29.  We  have 

ERF  i  (x)  =  erf(x),  ERFC  1  (0,  x)  =  erfc(x), 

where  the  subroutine  listings  for  these  functions  are  given  on  pages  F-12  to  F-15.  They  are  identical 
to  the  NSWC  (DL)  math  library  functions  ERF  and  ERFC  as  of  June  1980  which  are  based  on  the 
reference  below. 

Cody,  W.  J.,  Rational  Chebyshev  Approximations  for  the  Error  Function,  Mathematic*  of 
Computation,  v.  23  ( 1 969),  pp.  63 1-637. 
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SUBROUTINE  DZ  (  X,Y,N,ANS,A  , IOP  ) 
DIMENSION  X{1),Y<1),U(4),V(4) 

A=0. 

ANS=0. 

IF  (  N.NE.l  )  GO  TO  3031 
CALL  DZl  (  X,Y,N,ANS,IOP,A  ) 

RETURN 

3031  CONTINUE 

IF  (  N.LT.3  )  RETURN 
L=3 

U(1)=X(1) 

U  (2)  =X  (2) 

U  (3)  =X  (3) 

V(1)=Y(1) 

V{2) =Y (2) 

V  ( 3 )  =Y  { 3 ) 

3041  CONTINUE 

CALL  DZl  (U,V,3,ANSl,IOP,Al  ) 
A«A+A1 
ANS«ANS+ANS1 
3061  CONTINUE 
L^L+l 

IF  (  L.GT.N  )  RETURN 
U(2)=U(3) 

V(2)  ®V(3) 

U  {3}  »X  (L) 

V(3) «Y (L) 

GO  TO  3041 
END 
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SUBROUTINE  DZ-1 


Computes  P(Aj)  for  DZ 


CALL:  DZ-1  (x,  y,  N,  P,  IOP,  A)*, 
where: 

x  is  the  array  of  abscissas  of  the  points  specifying  a  simple  polygon  S. 

y  is  the  array  of  ordinates  of  the  points  specifying  a  simple  polygon  S.  x  and  y  are  dimen¬ 
sioned  at  N  +  1 . 

N  is  the  number  of  points  specifying  S. 

IOP  is  specified  by  the  user. 

IOP  =  1  for  3-decimal-digir  accuracy  for  P(S). 

IOP  =  2  for  6-decimal-digit  accuracy  for  PCS). 

IOP  =  3  for  9-decimal-digit  accuracy  for  P(S). 

P,  A  are  the  locations  where  the  values  of  P(S)  and  A(S)  are  returned,  respectively. 

For  each  angular  region  a  of  S  specified  by  R,  0lt  02 ,  DZ-1  computes  the  corresponding 
Drezner  arguments  m,  k.  p  as  indicated  in  (  2,  Eq.  60J.  Subroutine  PLAN  uses  Drezncr’s  algorithm 
to  determine  which  equation  of  (t]  is  used  to  find  P(a).  Functions  EQ  7,  EQ  8,  EQ  9  and  EQ  1 1  of 
PLAN  compute  P(a)  using  equations  7,  8,  9,  and  1 1,  respectively,  of  (t).  Subroutine  BPHI  uses 
equation  5  of  it  1  to  compute  P(a). 


*l)Z- 1  computes  P(S) by  Dreawi's  procedure  which  U  described  in  |2) . 

fZ.  Dt&uwi,  Compulsion  of  the  Biwiate  Normal  Intctmi,  Mathematics  of  Compulation,  v.  32(1 978),  pp,  277-279. 


G-5 


SUBROUTINE  DZl  {  X,Y,N, ANS, IOP, A  ) 
DIMENSION  X(l) ,Y(1) ,H(2) ,APH1(3) 
DATA  (  APH1(I) ,1=1,3  )  / 

1  2.02E-7 , 2.08E-13 ,2. 72E-19  / 

DATA  RT2  /  1.4142  13562  3731/ 

DATA  TWOPI/6 . 2831  85307  17958  / 

K  =  1 
ANS=0. 

IF  (  N.NE.l  )  GO  TO  3071 
W=X (2) -X (1) 

Z=Y (2) -Y (1) 

U=X (3) -X (1) 

V=Y (3) -Y (1) 

PSI1=(V*W-U*Z) 

IF  (  PSI1.GE.0.  )  GO  TO  3081 

ANS=+1. 

Tl=W 

W=U 

U=Tl 

T1=V 

V=Z 

Z*=T1 

GO  TO  3081 
3071  CONTINUE 

X (N+l) =X  (1) 

Y (N+l) =Y (1) 

CALL  SMP7  (  N, A,X,Y  ) 

IF  (  ABS(A) . LE.0.6512E-10  )  RETURN 
W=X (1) -X (N) 

Z=Y(1)-Y(N) 

U«X (2) -X (1) 

V=Y (2) -Y (1) 

3081  CONTINUE 

BGD1«SQRT(2. * (W*W+Z*2) ) 

BGD2*SQRT ( 2 . * (U*U+V*V) ) 

3151  CONTINUE 

B“.5*  (X  (K)  *X  (K)  +Y  (K)  *Y  (K) ) 

IF  (  B.GT. APH1 { IOP)  )  GO  TO  3155 

Tl*V*W-U*Z 

T2»U*W+V*Z 

PHIK=ATAN2  (Tl,T2  ) 

ANSlaPHIK/TWOPI 
GO  TO  3211 
3155  CONTINUE 

RTR= ( 2 . #  (W*V-U*Z))/(BGDl*BGD2) 

H(l)  ={-Y(K)*W+X(K)*Z)/BGDl 
H(2) =(-Y(K) *U+X(K) *V)/BGD2 
SGN=1. 
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IF  (  RTR.GE.O.  }  GO  TO  3161 

RTR=-RTR 

SGN=-1. 

T1=H{1) 

H  (1)  =H (2) 

H{2)=T1 
3161  CONTINUE 

AM  =-RT2*H (2) 

AK  =RT2*H (1) 

RHO= (-2. * (W*U+V*Z) ) / (BGD1*BGD2) 

IF  (  ABS(RHO) .LT. (1.-1.E-13)  )  GO  TO  3181 

IF  (  RHO.LT.O.  )  GO  TO  3171 

Tl=AM 

IF  (  AK.LE.AM  )  Tl=AK 
T2=-Tl/RT2 
ANS1=.5*ERFC1(0,T2  ) 

GO  TO  3191 
3171  CONTINUE 
ANS1=0 . 

IF  (  AK.LE.-AM  }  GO  TO  3191 

T1=-AK/RT2 

T2=AM/RT2 

ANSI*. 5* (ERFC1 (0,T.l) -ERFC1 (0 ,T2) ) 

GO  TO  3191 
3181  CONTINUE 

CALL  PLAN  (  AM  ,AK  ,RHO  , ANSI, IOP,RTR  ) 

3191  CONTINUE 

ANSl=SGN*ANSl 
3211  CONTINUE 

IF  {  K.NE.N  )  GO  TO  3651 
IF  (  N.NE.l  )  GO  TO  3631 
ANS« ABS ( ANS-ABS (ANSI)) 

RETURN 

3631  CONTINUE 

ANSeANS-ANSl 

IF  (  A.LT.O.  )  GO  TO  3641 
ANS=AN$+1. 

RETURN 

3641  CONTINUE 
ANS^ANS-l. 

RETURN 

3651  CONTINUE 
K=K+1 
KPl=K+l 
K=U 
Z=V 

U°X (KPl) -X (K) 

V-Y (KPl)-Y (K) 
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BGD1=BGD2 

BGD2=SQRT (2 . * (U*U+V*V) ) 

ANS=ANS-ANSl 

GO  TO  3151 

END 


SUBROUTINE  PLAN  (  H, AK,R,ANS, IOP  ,RTR  ) 
ANS=0 . 

IF  (  (H*AK*R) .GT.O .  )  GO  TO  3155 
IF  (  H.GT.O.  )  GO  TO  2031 
IF  {  AK.GT.O.  )  GO  TO  2021 
IF  (  R.GT.O.  )  GO  TO  2011 
ANS=BPHI(H,AK,R,IOP,RTR  ) 

RETURN 

2011  CONTINUE 

IF  (  AX.NE.O.  )  GO  TO  2061 
GO  TO  2023 
2021  CONTINUE 

IF  {  R.LT.O.  )  GO  TO  2041 
2023  CONTINUE 

ANS=EQ9(H,AK,R,IOP,RTR  ) 

RETURN 

2031  CONTINUE 

IF  (  AK.EQ.O.  )  GO  TO  2051 
2035  CONTINUE 

IF  (  AK.LT.Q.  )  GO  TO  2061 
2041  CONTINUE 

ANS=EQ7(H,AK,R,IOP,RTR  ) 

RETURN 

2051  CONTINUE 

IF  (  R.GT.O.  )  GO  TO  2061 
GO  TO  2041 
2061  CONTINUE 

ANS«EQ8(H,AK,R,IOP,RTR  ) 

RETURN 

3155  CONTINUE 

ANSbEQ11 (H, AK,R, IOP,RTR  ) 

RETURN 

END 


FUNCTION  EQ7  (H,AK,R,IOP,RTR  ) 
DATA  RT2/1.4142  13562  3731/ 
T=-H/RT2 
T1=-AK/RT2 

EQ7=BPHI(-H,-AK,R,I0P,RTR  ) 

1  +.5*  (ERFC1 (0 ,T)  +ERFC1  ( t;  , Tl) )  -1. 

RETURN 
END 


FUNCTION  EQ8  (H,AK, R, IOP, RTR  ) 

DATA  RT2/1.4142  13562  3731/ 

TI=-AK/RT2 

EQ8=-BPHI (-H, AK ,-R, IOP, RTR) +. 5*ERFC1 ( 0,T  ) 

RETURN 

END 


FUNCTION  EQ9  {H, AK,R,lOP,RTR  ) 

DATA  RT2/1.4142  13562  3731/ 

T<*-R/R’f2 

EQ9““BPRI (H,-AK,-R, I0P,RTR) +.5*ERFCl{0  * 

RETURN 

END 
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FUNCTION  EQU(K,AK,R,IOP,RTR  ) 

DATA  RT2/1.4142  13562  3731/ 

CST=SQRT (H*H-2. *R*H*AK+AK*AK  ) 

Tl=R*H-AK 

Cl=l. 

T2=SIGN{C1,H) 

T1=(T1*T2)  /CST 
T4=l. 

T3=H*AK 

T5=SIGN(T4,T3) 

TDEL= ( 1 . -  T5) *.25 

T3=R*AK-H 

Cl=l. 

T2=SIGN {Cl, AK) 

T3={T3*T2)  /CST 
RTR1= (RTR*ABS (H) ) /CST 
RTR3= (RTR*ABS ( AK) ) /CST 
IF  (  H.GT.O.  )  GO  TO  2031 
IF  {  Tl.GT.O.  )  GO  TO  2023 
T4=BPHI (H,0. ,T1, IOP,RTRl) 

GO  TO  2051 
2023  CONTINUE 

T4=EQ9(H,0.,Tl,IOP,RTRl) 

GO  TO  2051 
2031  CONTINUE 

IF  (  Tl.LT.Q.  )  GO  TO  2041 
T4«.5-BPHI(-H,0.,-T1,IOP,RTR1  ) 

GO  TO  2051 
2041  CONTINUE 
Cl=-H/RT2 

T4  *BPHI (-H«0. ,T1, IOP, RTRl ) - . 5  *ERF1 (C 1 ) 
2051  CONTINUE 

IF  (  AK.GT.0.  )  GO  TO  3031 
IF  (  T3.GT.0.  )  GO  TO  3023 
T6  *BPHI ( AK , 0 . , T3 , IOP , RTR3 ) 

GO  TO  3051 
3023  CONTINUE 

T6-EQ9(AK,0.,T3,IOP,RTR3  ) 

GO  TO  3051 
3031  CONTINUE 

IF  (  T3.LT.0.  )  GO  TO  3041 
T6*.S-BPHI(-AK,0.,-T3,IOP,RTR3  ) 

GO  TO  3051 
3041  CONTINUE 

C1“*-AK/RT2 

T6 =BPHI ( -AK , 0 . , T3 , IOP, RTR3  ) - * 5*ERF1 (Cl) 
3051  CONTINUE 

EQll*T4+T6-TDEL 
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RETURN 

END 


G-12 


1 

2 

3 

4 

1 

2 

3 

4 

1 

2 

3 

4 

1 

2 

3 

4 


1 


3351 

3361 


FUNCTION  BPHI  (  H,AK,R  ,I0P,RTR  ) 
DIMENSION  A(21) ,X (21) ,LLO(6) ,LHI (6) 
DIMENSION  EPS 1(11  ) 


DIMENSION  EPS3 (11) 

DATA  (  A( I) ,1=1,8  )  / 
4.4602  97704  66658E-1, 
4.3728  88798  77644E-2, 
3.9233  10666  52399E-1, 
3.3246  66035  13439E-2, 
DATA  (  X(I) ,1=1,8  )  / 
1.9055  41497  98192E-1, 
1.7997  76578  41573E+0, 
4.8281  39660  46201E-1, 

1.7797  29418  52026E+0, 


3.9646  82669 
2.4840  61520 
2.1141  81930 
8.2485  33445 


98335E-1, 
28443E-1, 
76057E-1, 
15628E-4  / 


8.4825  18675  44577E-1, 
1.0024  21519  68216E-1, 
1.0609  49821  52572E+0, 
2.6697  60356  08766E+0  / 


DATA  (  A(I) , 1=9 ,16  )  / 

1.3410  91884  53360E-1,  2.6833  07544  72640E-1, 
2.7595  33979  88422E-1,  1.5744  82826  13790E-1, 
4.4814  10991  74625E-2 ,  5.3679  35756  02526E-3, 
2.0206  36491  32407E-4,  1.1925  96926  59532E-6  / 
DATA  (  X(I) ,1=9,16  )  / 

5. 2978  64393  18514E-2,  2,6739  83721  67767E-1, 
6.1630  28841  82402  E-l,  1.0642  46312  11623E+0, 
1.5888  55862  27006E+0,  2.1839  21153  09586E+0, 
2.8631  33883  70808E+0 ,  3.6860  07162  72440E+0  / 
DATA  (  EPSl(I) ,1=1,3  )  /  -8. ,-12. ,-20.  / 

DATA  PI  /  3.1415  92653  58979  / 

DATA  (  LLO (I) ,1=1,3  )  /  1,4,9  / 

DATA  (  LHKI), 1=1, 3  )  /  3,8,16  / 

DATA  RT2  /  1.4142  13562  3731/ 

DATA  (  EPS3(I) ,1=1,3  )  /  2.E-5,2.E-7,2.E-10  / 
ILO=LLO ( IOP) 

IHI=LHI ( IOP) 

EPS=£PS1 ( IOP) 

CST=RT2*RTR 


BPHI=0. 

H1=H/CST 
AK1=AK/CST 
SUM=0 . 

DO  3361  I=ILO, IHI 
SUM1=0. 

DO  3351  J=ILO, IHI 

Tl=Hl* (2 . *X { I) -HI) +AK1* (2. *X ( J) -AKl) 
+2.*R*(X(I)-H1)*(X(J)-AK1) 

IF  (  T1.LT.EPS  )  GO  TO  3351 
S'JMI=SUM1+EXP (Tl)  *A(J) 

CONTINUE 

SUM=SUM4-A( I)  *SUM1 
CONTINUE 
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BPHI=(SUM*RTR)/PI 

RETURN 

END 
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